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Abstract — In this paper we present a generic framework for 
the asymptotic performance analysis of subspace-based param- 
eter estimation schemes. It is based on earlier results on an 
explicit first-order expansion of the estimation error in the signal 
subspace obtained via an SVD of the noisy observation matrix. 
We extend these results in a number of aspects. Firstly, we 
demonstrate that an explicit first-order expansion of the Higher- 
Order SVD (HOSVD)-based subspace estimate can be derived. 
Secondly, we show how to obtain explicit first-order expansions 
of the estimation error of arbitrary ESPRIT-type algorithms and 
provide the expressions for R-D matrix-based and tensor-based 
Standard ESPRIT as well as Unitary ESPRIT. Thirdly, we derive 
closed-form expressions for the mean square error (MSE) and 
show that they only depend on the second-order moments of the 
noise. Hence, we only need the noise to be zero mean and possess 
finite second order moments. Additional assumptions such as 
Gaussianity or circular symmetry are not needed. Fourthly, we 
investigate the effect of using Structured Least Squares (SLS) to 
solve the overdetermined shift invariance equations in ESPRIT 
and provide an explicit first-order expansion as well as a closed- 
form MSE expression. Finally, we simplify the MSE for the 
special case of a single source and compute the asymptotic 
efficiency of the investigated ESPRIT-type algorithms in compact 
closed-form expressions which only depend on the array size and 
the effective SNR. 

Our results are more general than existing results on the 
performance analysis of ESPRIT-type algorithms since (a) we do 
not need any assumptions about the noise except for the mean to 
be zero and the second-order moments to be finite (in contrast 
to earlier results that require Gaussianity and/or second-order 
circular symmetry); (b) our results are asymptotic in the effective 
SNR, i.e., we do not require the number of samples to be large 
(in fact we can analyze even the single-snapshot case); (c) we 
present a framework that incorporates the SVD-based and the 
HOSVD-based subspace estimates as well as Structured Least 
Squares in one unified manner. 
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I. Introduction 

IGH resolution parameter estimation from R- 
dimensional (R-D) signals is a task required for a variety 
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of applications, such as estimating the multi-dimensional 
parameters of the dominant multipath components from 
MIMO channel measurements lfl2l . which may be used for 
geometry-based channel modeling. Other applications include 
radar If24l . wireless communications |f20l , sonar, seismology, 
and medical imaging. In ifTTl . we have shown that in the R-D 
case (R > 2), tensors can be used to store and manipulate the 
R-D signals in their native multidimensional form. Based on 
this idea, we have proposed an enhanced tensor-based signal 
subspace estimate as well as ESPRIT-type algorithms based 
on tensors in IfTTl . Their superior performance was shown 
based on Monte-Carlo simulations. 

In this paper we present a framework for the analytical 
performance assessment of subspace-based parameter estima- 
tion schemes and apply it to derive a first-order perturbation 
expansion for the tensor-based subspace estimate. Moreover, 
we find first-order expansions for the estimation error of 
ESPRIT-type algorithms and derive generic mean square error 
(MSE) expressions which only depend on the second-order 
moments of the noise and hence do not require Gaussianity 
or circular symmetry. This approach allows to assess the gain 
from using tensors instead of matrices analytically in order to 
determine in which scenarios it is particularly pronounced. We 
apply the framework for the analysis of R-D Standard ESPRIT, 
R-D Unitary ESPRIT, R-D Standard Tensor-ESPRIT and R-D 
Unitary Tensor-ESPRIT. Moreover, we investigate the effect 
of using Structured Least Squares (SLS) for the solution of 
the invariance equations. Finally, we present simplified MSE 
expressions for the special case of a single source impinging 
on a Uniform Linear Array (ULA) as well as a Uniform Rect- 
angular Array (URA) and observed under circularly symmetric 
white noise. These expressions only depend on the effective 
Signal to Noise Ratio (SNR) and the array size and the allow 
to compute the asymptotic efficiency in closed-form. 

Analytical performance assessment of subspace-based pa- 
rameter estimation schemes has a long standing history in 
signal processing. Shortly after the publication of the most 
prominent candidates, MUSIC 03 and ESPRIT OJJ, ana- 
lytical results on their performance have appeared. The most 
frequently cited papers are fBl for the MUSIC algorithm 
and [28 1 for ESPRIT. However, many follow-up papers exist 
which extend the original results, e.g., J26l, Q, (23, EH, 
(23 > a nd many others. However, these results have in common 
that they all go more or less directly back to a result on the 
distribution of the eigenvectors of a sample covariance matrix 
from (3. 
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In contrast to these results, in IfTTl an entirely different 
approach was proposed, which provides an explicit first-order 
expansion of the subspace of a desired signal component 
if observed superimposed by a small additive perturbation. 
This approach has a number of advantages compared to A3- 
Firstly, |3 is asymptotic in the sample size TV, i.e., the result 
becomes only accurate as the number of snapshots TV is very 
large, whereas IfTTl is asymptotic in the effective SNR, i.e., it 
can be used even for TV = 1 as long as the noise variance is 
sufficiently small. Secondly, (3 requires strong Gaussianity 
assumptions, not only on the perturbation (i.e., the noise), 
but also on the source symbols. Since IfTTl is explicit, no 
assumptions about the statistics of either the desired signal 
or the perturbation are needed. Note that it has recently been 
shown that 13 can be extended to the non-Gaussian and 
the non-circular case in (5). However, the large sample size 
assumption is still needed. Thirdly, the covariance expressions 
from (3 are much less intuitive than the expansion from ifTD 
which shows directly how much of the noise subspace "leaks 
into" the signal subspace due to the erroneous estimate. 
Finally, the expressions involved in 13 are quite complex and 
tough to handle, whereas IfTTl requires only a few terms which 
appear directly as block matrices of the SVD of the noise-free 
observation matrix. 

Due to these advantages we clearly favor IfTTl as a starting 
point. The authors in IfTTl have already shown that their results 
on the perturbation of the subspace can be used to find a first 
order expansion for the MUSIC, the Root-MUSIC, the Min- 
Norm, the State-Space-Realization, and even the ESPRIT algo- 
rithm. However, they only considered 1-D Standard ESPRIT. 
We extend their work by considering multiple dimensions 
(R-D ESPRIT), by incorporating forward-backward-averaging 
(for Unitary ESPRIT), by considering the tensor-based sub- 
space estimate (for Standard and Unitary Tensor-ESPRIT), 
by investigating the effect of using Structured Least Squares 
(SLS) to solve the invariance equation instead of the Least 
Squares (LS) solution used in IfTTl . and by providing generic 
mean square error (MSE) expressions of the resulting esti- 
mation errors in these cases. Note that our MSE expressions 
depend on the second-order moments of the noise only. Hence 
we only assume it to be zero mean (due to the asymptotic 
nature of our performance analysis), but do not require it 
to be Gaussian distributed, white, or circularly symmetric. 
This is a particularly attractive feature of our approach with 
respect to different types of preprocessing which alters the 
noise statistics, e.g., spatial smoothing (which yields spatially 
correlated noise) or forward-backward averaging (which an- 
nihilates the circular symmetry of the noise). Since we do 
not require spatial whiteness or circular symmetry, our MSE 
expressions are directly applicable to a wide range of ESPRIT- 
type algorithms. 

There have been other follow-up papers based on IfTTl . 
For instance, l40l provides a first-order and second-order 
perturbation expansion which can be seen as a generalization 
of IfTTl . In lfT9l the authors show that there is also a first- 
order contribution of the perturbation of the signal subspace 
which lies in the signal subspace (which |40) and IfTTl have 
argued to be of second order and hence negligible). Note that a 



perturbation expansions for the signal subspace and null space 
projectors based on the sample covariance matrix is provided 
in lfl5l where the expansion up to an arbitrary order is derived 
based on a recurrence relation. A mean square error expression 
for Standard ESPRIT is provided in lfT"8l . however, it does not 
generalize easily to the tensor case and it assumes circular 
symmetry of the noise. Note that the latter assumption implies 
that it is not applicable to Unitary ESPRIT, since Forward- 
Backward Averaging annihilates the circular symmetry of the 
noise. Moreover, other authors have studied the asymptotical 
performance of ESPRIT, e.g., J35), J6) where harmonic re- 
trieval from time series is investigated and MSE expressions 
for a large number of snapshots as well as MSE expressions for 
a high SNR are derived. Note that we find MSE expressions 
compatible to (61 by only assuming a high effective SNR, 
i.e., either the number of snapshots or the SNR can tend to 
infinity. Interestingly, J35), J6l also consider the special case 
for a single source. However, the expressions provided there 
are specific to harmonic retrieval from time series and they are 
not compared to the corresponding Cramer-Rao Bound. Some 
analytical results on the asymptotic efficiency of MUSIC, 
Root-MUSIC, ESPRIT, and TLS-ESPRIT are, among others, 
presented in (27), l28l . l29l . and ll25l . respectively. However, 
these results are asymptotic in the number of snapshots TV and 
sometimes even in the number of sensors M. The asymptotic 
equivalence of LS-ESPRIT, TLS-ESPRIT, Pro-ESPRIT, and 
the Matrix Pencil method has also been shown, see for instance 
|fl3l . Overall, in the matrix case, the number of existing results 
is quite large, since the underlying methods have been known 
for more than two decades. However, concerning the tensor 
case and the incorporation of Structured Least Squares, the 
existing results are much more scarce. 

In the tensor case a first-order expansion for the HOSVD 
has been proposed in [4|. However, it is not suitable for 
our application since it does not consider the HOSVD-based 
subspace estimate but the subspaces of the separate n-mode 
unfoldings and their singular values. Moreover, the pertur- 
bation is modeled via a single scalar real-valued parameter 
e. A first-order expansion for the best rank-(i?i, R2, R3)- 
expansion is provided in 0. However, again, it is not di- 
rectly applicable for analyzing the HOSVD-based subspace 
estimate as it investigates the approximation error of the entire 
tensor. Consequently, our approach to analyze the HOSVD- 
based subspace estimate based on the link to the SVD-based 
subspace estimate via a structured projection is entirely novel. 
Moreover, the application of these results to find the analytical 
performance of Tensor-ESPRIT-type algorithms is novel as 
well. It is a particular strength of the framework we use 
that many extensions and modifications of ESPRIT are easily 
incorporated, e.g., Forward-Backward-Averaging or Structured 
Least Squares. 

This paper is organized as follows: The notation and the data 
model are introduced in Sections [TT] and [HI] respectively. The 
subsequent Section [TV] reviews the first-order perturbation of 
the matrix-based subspace estimate and presents the extension 
to the tensor case. The performance analysis of ESPRIT-type 
algorithms is shown in [VJ Numerical results are presented in 
Section [VTl before drawing the conclusions in Section fVHl 
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II. Notation 

In order to facilitate the distinction between scalars, matri- 
ces, and tensors, the following notation is used: Scalars are 
denoted as italic letters (a, b, . . . , A, B, . . . , a, j3, . . .), column 
vectors as lower-case bold-face letters (a, b, . . .), matrices 
as bold-face capitals (A, B, . . .), and tensors are written as 
bold-face calligraphic letters (A, B, . . .)■ Lower-order parts 
are consistently named: the -element of the matrix A, 
is denoted as Ojj and the (i, j, fc)-element of a third order 
tensor B as bij^k- 

We use the superscripts T , H ,* ~ 1 , + for transposition, Her- 
mitian transposition, complex conjugation, matrix inversion, 
and the Moore-Penrose pseudo inverse of a matrix, respec- 
tively. The trace of a matrix A is written as Tr (A). Moreover, 
the Kronecker product of two matrices A and B is denoted as 
A g) B and the Khatri-Rao product (column-wise Kronecker 
product) as AoB. The operator vcc {A} stacks the column of 
a matrix A £ C Mx N into a column vector of length M- N X 1. 
It satisfies the following property 

vec{A-X-B} = (B T <g> A) -vec{X}. (1) 

An n-mode vector of an (Ii x I2 X . . . x In) -dimensional 
tensor A is an /„ -dimensional vector obtained from A by 
varying the index i n and keeping the other indices fixed. 
Moreover, a matrix unfolding of the tensor A. along the n-fh 
mode is denoted by [A] ^ and can be understood as a matrix 
containing all the n-mode vectors of the tensor A. The order 
of the columns is chosen in accordance with ||4). 

The outer product of the tensors A £ c 7lX/2X - x/jv and 
B e C JlXj2X --- XJju is given by 

C = Ao B e C /iX --- x7jvxJiX --- xJm , where 

7 

Ci lt i 2 ,...,i N ,ji,32,— ,jM °»i,»2 *w ' 31 ,32 ,■ ■ ■ ,3 m ' 

In other words, the tensor C contains all possible combinations 
of pairwise products between the elements of A and B. This 
operator is very closely related to the Kronecker product 
defined for matrices. 

The n-mode product of a tensor A £ Qiixi 2 x...xi N an( j 
a matrix U £ C' , ™ xJ " along the n-th mode is denoted as 
B = A x „ U and defined via 

B = Ax n U &[B] {n) = U-[A\ n) , (3) 

i.e., it may be visualized by multiplying all n-mode vectors 
of A from the left-hand side by the matrix U. Note that the 
n-mode product satisifies 

(A X r U r ) Xr V r = A X r (V r ■ U r ) (4) 

[AxiUx... x R U R ] {r) = U r ■ [A] {r) ■ 
(U r+1 (8 . . . <8 U R (8 Ut ® . . . ® Lf r _i) T (5) 

for A£ C hxhx - xlN , U r £ C J - XI - and V r £ C K - Xj -. 
The higher-order SVD (HOSVD) § of a tensor A £ 

c / 1 x/ 2 x...x/„ j s given by 

A = S xi U\ x 2 C7 2 . . . xjv Un, (6) 

where S £ C JlXj2X '" xiN is the core tensor which satisfies 
the all-orthogonality conditions JH and U n £ C InXl ™, n = 



1, 2, . . . , N, are the unitary matrices of n-mode singular vec- 
tors. 

We also define the concatenation of two tensors along the 
n-th mode via the operator [A i_i n B] . The Euclidean (vector) 
norm, the Frobenius (matrix) norm, and the Higher-Order 
Frobenius (tensor) norm are denoted by ||a|| 2 , ||A|| F , and 
|*4|| H , respectively. All three norms are computed by taking 
the square-root of the sum of the squared magnitude of all the 
elements in their arguments. 

The matrix K^i x n denotes the commutation matrices ||2T1| 
which satisfy 

K M xn ■ vec {a t | = vcc {A} (7) 

Kj IxP -{A®B)-K NxQ = B(g)A (8) 

for A £ C MxN , B £ C PxQ . 

A projection matrix onto the column space of a matrix 
A £ C Mxr is denoted as Pr A = A A + £ C MxM and 
its orthogonal complement by FV^ = 1 14 — Pr^. Note that 
for r = 1, i.e., A = a, this matrix can be computed as 

I V _ a a H 
rL a — a H. a - 

A p x p matrix Q p is called left-FI-real if TJ p ■ Q* = Q p , 
where Il p is the p x p exchange matrix with ones on its 
antidiagonal and zeros elsewhere. The special set of unitary 
sparse left-II-real matrices introduced in J5) is denoted as 
Qp. Furthermore, a matrix X £ <C MxN is called centro- 
Hermitian if II m ■ X* ■ 11^ = X. The vector denotes the 
k-th column of an identity matrix. 

III. Data model 

A. Matrix-based and tensor-based data model 

The observations are modeled as a superposition of d 
undamped exponentials sampled on an i?-dimensional grid of 
size Mi x Mi x ... x Mr at N subsequent time instants iflOl . 
The measurement samples are given by 

d R 

x t - Vs-fi ) TTe J ' (m '- 1) -^ > +n 

i=l r=l 

(9) 

where m r = 1,2,..., M r , n = 1, 2, . . . , N, Si(t n ) denotes the 
complex amplitude of the i-th exponential at time instant t n , 

(r) 

fi\ symbolizes the spatial frequency of the i-th exponential 
in the r-th mode for i = 1,2, ... ,d and r = 1, 2, . . . , R, 
and n mi)m2i ... jTOfl) t n represents the zero mean additive noise 
component inherent in the measurement procesfl In the 
context of array signal processing, each of the i?-dimensional 
exponentials represents one planar wavefront and the complex 
amplitudes Si(t n ) are the symbols. It is our goal to estimate 

(r) 

the spatial frequencies ^\ for r = 1, 2, . . . , R, i = 1, 2, . . . , d 
and their correct pairing. 

In order to arrive at a more compressed formulation of the 
data model in (O we collect the samples x mi ,m 2 ,...,m R ,t n into 
one array. As our signal is referenced by R + 1 indices, the 

1 Note that equation j9) assumes a uniform sampling in the spatial domain. 
However, this assumption can be relaxed to more generic geometries as long 
as they feature shift invariances and can be constructed as the outer product 
of R one-dimensional sampling grids. 
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most natural way of formulating the model is to employ an 
(R+ l)-way array X G C MlXM2 ~ xM * xJV which contains 
x mumi2l ...,mR,t n for m r = 1, 2, . . . , M r , r = 1, 2, . . . , R, and 
n = 1, 2, . . . , N, We can then conveniently express X as ifTTl 

X = Ax R+1 S T +M, (10) 

where S G <C dxN contains the amplitudes Sj[n] and 
A/" G C MlXAr2 - xMRXJV collects all the noise samples 
n m 1 ,m 2 ....,m R ,t n m the same manner as X. Finally, A G 
pMixJ/ 2 ...xA/ H x(i j s re f erre( j as the "array steering tensor" 

ffTTl . It can be expressed by virtue of the concatenation 
operator via 

<A= [^4-1 1 — lit-j-i -4-2 1 — ir+1 ■ ■ ■ i—iR+i<Ad] (11) 

A = aU(^) o a^(^) 0...0 a^(^) G C M ^ M **~ 



e C PlX 

Af r xp r 



R+l 



•<Nxd 



where S 

- [s] 

and [/,, G 

are the matrices of dominant r-mode singular vectors. Here, 
rank I [Xq] , s ) represents the r-rank of the noise-free 



..xp R xd j s ^ truncated, core tensor 
for r = 1,2, ...,R, U 



Pr 



observation tensor Xo = Ax r + i S . Based on (TOI l. a tensor- 
based subspace estimate can be defined as 



U = S x 1 U 1 



TT*-' 

XR Uj 



XR+l^s ■ 



(16) 



t: 1 



IS 



where (/i 



■>M r xl 



represents the array steering vector 



of the i-th source in the r-th mode. 

An alternative expression for the array steering tensor is 
given by 

A = X R+1 , d x 1 x 2 A {2) ... x R A iR \ (12) 



Note that the (R + l)-mode multiplication with 
introduced in addition to its original definition in IfTTl since 
it simplifies the notation we need at this point and it has no 
impact on the subspace estimation accuracy. Here S s refers to 
the diagonal matrix containing the d dominant singular values 
on its main diagonal. 

~[s] *[s] . 

Note that U satisfies U ~ A x r + i T for a non-singular 
matrix T € C dxd . Based on Zi' , an improved signal subspace 



aM(/4 r) ) 



a 



G C M - Xd is 



where A^ 

referred to as the array steering matrix in the r-th mode. 

The strength of the data model in ( TfOl i is that it represents 
the signal in its natural multidimensional structure by virtue 
of the measurement tensor X. Before tensor calculus was 
used in this area, a matrix-based formulation of ( TTOb was 
needed. This requires stacking some of the dimensions into 
rows or columns. A meaningful definition of a measurement 
matrix X is to apply stacking to all "spatial" dimensions 
1,2, ...,R along the rows and align the snapshots n = 
1,2, ...,N as the columns. Mathematically, we can write 
X = [X]J R+1) G C MxN , where M = l\f=i M r- Applying 
this stacking operation to (flOb . we arrive at the matrix-based 
data model IfTOl 



X = A-S + N. 



Here, A = [A] 



(R+l) 



•>A/xd 



and N = [Af] 



(13) 

±) G C MxN . 



Note that A is highly structured since it satisfies 

A=[A]J R+1) =A^oA ( Vo...oA^. 



(14) 



B. Subspace estimation 

The first step in all subspace-based parameter estimation 
schemes is the estimation of a basis for the signal subspace 
from the noisy observations. In the matrix case, this can, for 
instance, be achieved by a truncated SVD of X. Let XJ S G 
C Mxd be the matrix containing the d dominant left singular 
vectors of X. Then the column space of XJ S is an estimate 
for the signal subspace spanned by the columns of A and we 
can write A « U s ■ T for a non-singular matrix T G <C dxd . 

A tensor-based extension of this subspace estimate was 
proposed in IfTTl . To this end, let the truncated HOSVD of 
X be given by 



X rj S x 1 U 1 



x R U r x r+1 U 



R+l J 



(15) 



estimate is given by the matrix 



U 



(R+l) 



G C 



Mxd 



IV. Perturbations of the subspace estimates 

A. Review of perturbation results for the SVD 

Let us first review the results from [17| which are relevant 
to the discussion in this section. Let Xq = A ■ S G C MxN 
be a matrix containing the noise-free observations such that 
X = Xq + N where N represents the undesired perturbation 
(noise). 

The SVD of Xo can be expressed as 



X = [U s U n ] 



where the columns of U s G 



S s 




[V s V n ] H , (17) 



provide an orthonormal 



basis for the signal subspace which we want to estimate. 
Moreover S s = diag ( \a\, a-2, ■ ■ ■ , Cd] ) G R dxd contains the 
d non-zero singular values on its main diagonal. We find an 
estimate for U s by computing an SVD of the noisy observation 
matrix X which can be expressed as 



x = [u s u n ] 











(18) 



where the "hat" denotes the estimated quantities. We can write 
U s = U s + AU S , where AU S represents the estimation error. 
At this point we are ready to state the main result on the first 
order perturbation expansion of AU S from lfT7l 



AU S = U n 

r 



C{A 2 }, where A = ||iV|| and 



U^ ■ N- V s ■ S7 1 G c (M - d)xd (19) 

Here |.| represents an arbitrary sub-multiplicativ^| norm, 
e.g., the Frobenius norm. Equation ( fT9l shows the first order 
expansion of the signal subspace estimation error AU S in 
terms of the noise subspace U n , i.e., how much of the noise 
subspace "leaks into" the signal subspace due to the estimation 
errors from the perturbation TV. Since it is explicit in N it 

2 A matrix norm is called submultiplicative if \\A ■ B\\ < • ||B|| for 
arbitrary matrices A and B. 
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makes no assumptions about the statistics of N, in fact, it is 
purely deterministic. 

The expansion < fT~9b only models the leakage of the noise 
subspace into the signal subspace. That is to say, the pertur- 
bation of the particular basis (the columns of U s ) is ignored. 
While for subspace-based parameter estimation schemes this 
is indeed sufficient since the particular choice of the basis is 
irrelevant, there are other applications where this term matters. 
For instance, in a communication system where the channel 
is decomposed into its individual eigenmodes and one or 
several of these eigenmodes are used for transmission, such 
errors have a major impact. Therefore, other authors have 
extended ( fT9l ) to take this error term into account. For instance, 
in |fl9l the authors provide the following expansion 



AU S 



= u n 

= £>G 



J7 s -r s + C{A 2 }, where 



(20) 



Uf ■ N ■ V s • E s + S s • V? • iV H • U s ) e C p 



Here, the matrix D is defined as 

k = i 



\D] 



(k,£) 



for k, 



1,2,.. 



(21) 



Equation (120b additionally shows the perturbation of the in- 
dividual singular vectors via the term U s ■ T s . This term can 
be dropped for the evaluation of subspace-based parameter 
estimation schemes since for these, the particular choice of the 
basis is irrelevant. Therefore we do not consider it in SectionM 
where ESPRIT-type algorithms are investigated. However, we 
show its impact in the simulation results in Section |VT] where 
the subspace estimation accuracy is evaluated. 

B. Extension to the HOSVD-based subspace estimate 

As we have shown in Section UlI-BI in the multidimensional 
case, an improved signal subspace estimate can be computed 
via the HOSVD of the measurement tensor X. Since the 
HOSVD is computed via SVDs of the unfoldings, we can 
apply the same framework to find a perturbation expansion of 
the HOSVD-based subspace estimate. In order to distinguish 
unperturbed from estimated (perturbed) quantities we express 
X as X = Xq + A^where Xq = A. Xr + i S T is the unper- 
turbed observation tensor. The SVD of the 7--mode unfoldings 
of X and Xq are then given by 



Xo] (r) =[uW UW] 



[X] 



(r) 



uv u: 



diag 





r 

o 



[Vf V^f (22) 



si 



V 



(23) 



where S' sl 

1, 2, . . . , R. Note that since 
we can apply 

At/[ s] = U 

ri n i = u 



>) Jr) 



,o~2 ')•••) o"^] J an d r = 
and ([23j are in fact SVDs, 

= U [ * ] + A[/[ s] where 

rj s] +e>{A 2 }, (24) 



and find U r 

■ r[. nl + ul 



(r) 



D r 



VI 



SN . • Mf r) • C/N) 







fork,£ = 1,2,.. 



Our goal is to use the perturbation of the r-mode unfoldings 
to find a corresponding expansion for the HOSVD-based sub- 
space estimate introduced in Section IIII-BI This is facilitated 
by the following theorem: 

Theorem 1. The HOSVD-based subspace estimate 

T 



u 



(R+l) 



defined in ( 1161 l is linked to the SVD-based 



subspace estimate U s via the following algebraic relation 



u 



where T r £ 



(R+l) 



•>M r x M. 



■ U s , (25) 



r represent estimates of the projection 



matrices onto the r-spaces of Xq, which are computed via 

Proof: Relation (EU was shown in [30) for R = 2. The 
proof for R > 2 proceeds in an analogous manner and is 
presented in Appendix lAl 

Equation (l25T > shows that a perturbation expansion for 

-[s]l T 

U can be developed based on the subspaces of all 

J (R+l) 

R+l unfoldings, as the core tensor is not needed for its 
computation. The result is shown in the following theorem: 

Theorem 2. The HOSVD-based signal subspace estimate can 

r.r.nT r . r„n T 



be written as 



U 



AU 



T 

(R+l) 



= U S 

(R+i) 
(Ti ® T 2 (8 . 



AU 



(R+i) 



where 



AU™ ■ Ul 



T R ) ■ AU S 

5 T 2 ® ■ ■ ■ (8 T R 

5 . . . ® Tr 



0{A 2 }, 



AC/| ] • uf 



■U s 

■u s 
u s 

(26) 



the SVD-based signal subspace perturbation AU S is given 
by ( 120b and the perturbation of the r-space can be computed 
via 



AC/M = C/M • if! - C/W • (7[, n]H • [A0 (r) • VI s1 • EM \ 

(27) 



Proof: cf. Appendix IB"! 

Note that while AU S in general contains both perturbation 
terms U n ■ T n and U s ■ T s , for AU [ * ] the term U^ ] ■ r|f ] 
cancels. This is not surprising since the r-mode subspaces 
enter d25l l only via projection matrices for which the choice 
of the particular basis is irrelevant. 
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V. Asymptotical analysis of the parameter 

ESTIMATION ACCURACY 

A. Review of perturbation results for the 1-D Standard ES- 
PRIT 

In |[T7| the authors point out that once a first order expansion 
of the subspace estimation error is available it can be used to 
find a corresponding first order expansion of the estimation 
error of a suitable parameter estimation scheme. One of the 
examples the authors show is the 1-D Standard ESPRIT 
algorithm using LS, which we use as a starting point to 
discuss various ESPRIT-type algorithms in this section. In the 
noise-free case, the shift invariance equation for 1-D Standard 
ESPRIT can be expressed as 



Ji • U s ■ * = J 2 • U s 



(28) 



where J \,Ji £ 



pM (Bel> xM 



are the selection matrices that 



select the M t > scl '> elements from the M antenna elements 
which correspond to the first and the second subarray, re- 
spectively. Moreover, ^ = Q ■ • Q , where 3? = 
diag ( [e^ 1 , e JtJld ]) e £ dxd contains the spatial fre- 

quencies /ifc, k = 1, 2, . . . , d that we want to estimate. There- 
fore, fik = axg (EVfc i.e., the fc-th spatial frequency is 
obtained from the phase of the fc-th eigenvalue (EVfc {•}) of 
*. 

In presence of noise, we only have an estimate U s of the 
signal subspace U s . Consequently, (EHT l does in general not 
have an exact solution anymore. A simple way of finding 
an approximate if? is given by the LS solution which can be 
expressed as 



LS 



Ji-U t 



■J 2 -U s 



(29) 



To simplify the notation we skip the index "LS" for the 
remainder of this section (since only LS is considered) and 
pick it up again in the next section where we expand the 
discussion to SLS. 

For the estimation error of the fc-th spatial frequency cor- 
responding to the LS solution from ( f29l . ifTTl provides the 
following expansion 

A^ k =Im{p£ ■ {J 1 ■ U s )+ ■ [J 2 /X k - J x ] 

■AU s -q k }+(D{A 2 } (30) 

where X k = c mk and q k is the fc-th column of Q. Moreover, 

represents the fc-th row vector of the matrix P — 
Note that AU S can be expanded in terms of the perturbation 
term TV by directly using the expansion (fT9l . The additional 
term from (l20l i is not needed as it is irrelevant for the 
performance of ESPRIT. 

B. Extension to R-D Standard (Tensor- )ESPRIT 

The previous result from IfTTl on the first order perturba- 
tion expansion of 1-D Standard ESPRIT using LS is easily 
generalized to the R-D case. The reason is that for R- 
D LS-based ESPRIT, the R shift invariance equations are 



solved independently from each otheJH Hence, the arguments 
from ifPTl are readily applied to all modes individually and we 
directly obtain a first order expansion for the estimation error 
of the fc-th spatial frequency in the r-th mode 



A/4 r) =Im{ 
• AU 



T 
Pk 



fir) 
J 1 



U K 



(r) 
k 



fir) 
J 1 



Qk 



where J 1 , 



} +0{A 2 } 



M. M (^l) xM 



(31) 



are the effective R- 



D selection matrix for the first and the second subarray 
in the r-th mode, respectively. They can be expressed as 



I n r n -J 1 M n ® J 



Ar) 



1,2, ... ,R, where J\ ' G 

r(sel) 



for I = 1,2 and 

are the selection 
elements belonging to the 



M (=el) xM 



matrices which select the Mr 
first and the second subarray in the r-th mode, respectively. 

Since this expansion for R-D Standard ESPRIT is explicit 
in the perturbation of the subspace estimate and R-D Stan- 
dard Tensor-ESPRIT only differs in the fact that it uses the 
enhanced HOSVD-based subspace estimate, we immediately 
conclude that a first order perturbation expansion for R-D 
Standard Tensor-ESPRIT is given by 



A/4 r) =Im 



{pI • {jV 



u, 



f(r) n (r) _ fir) 



AU 



■q k )+0{A 2 }. 



(32) 



An explicit expansion of A(j£' in terms of the noise tensor 
Af is obtained by inserting the previous result d26l ). 

C. Mean Square Errors 

As it has been said above, the advantage of the first order 
perturbation expansion we have discussed so far is that it 
is explicit in the perturbation term N and hence makes 
no assumptions about its distribution. However, it is often 
also desirable to know the mean square error if a specific 
distribution is assumed and the ensemble average over all 
possible noise realizations is computed. 

In the sequel we show that the mean square error only 
depends on the second-order moments of the noise samples. 
Hence, we can derive the MSE as a function of the covariance 
matrix and the pseudo-covariance matrix only assuming the 
noise to be zero mean. We neither need to assume Gaus- 
sianity nor circular symmetry. For simplicity we consider the 
special case R = 2 for Standard Tensor-ESPRIT, however, 
a generalization to a larger number of dimensions is quite 
straightforward. 

3 If the shift invariance equations are solved completely independently, the 
correct pairing of the parameters across dimensions has to be found in a 
subsequent step. This is often avoided by computing the LS solutions for 
>f' r ) independently but then performing a joint eigendecomposition of all 
R dimensions to yield 4>' r '. This step is not included in the performance 
analysis presented in this section, since no performance results on joint 
eigendecompositions are available and it appears to be a very difficult task. 
Moreover, this step has indeed no impact on the asymptotic estimation error 
of the spatial frequencies for high SNRs since the eigenvectors become 
asymptotically equal. As shown in 1161 , the impact of the perturbation of 
the eigenvectors is of second-order and can hence be ignored in a first-order 
perturbation analysis. 
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Theorem 3. Assume that the entries of the perturbation term 
N or Af are zero mean random variables with finite second- 
order moments described by the covariance matrix R nn = 
E {n ■ n H } and the complementary covariance matrix C nn = 

E{n-n T } for n = vcc{7V} = vcc | [A/]^ |. Then, the 
first-order approximation of the mean square estimation error 
for the k-th spatial frequency in the r-th mode is given by 



E 



Rc 



f (r) 



W* mat ■ R 



T 
nn 



w 



be shown that this step has no impact on the performance for 
high SNRs. Therefore, the asymptotic performance of Unitary- 
ESPRIT-type algorithms is found once Forward-Backward- 
Averaging is taken into account. 

Forward-Backward-Averaging augments the N observations 
of the sampled R-D harmonics by N new "virtual" observa- 
tions which are a conjugated and row- as well as column- 
flipped version of the original ones [9|. This can be expressed 
in matrix form as 



w 



for R-D Standard ESPRIT and 

2~ 



2 \ r k 



W 



■ n T 

ten -^iin 



o{Tr(R nn f} 
(33) 



(fba 



= x n 



M 



x* -n N ] e c 



Mx2N 



Inserting X = Xq 

X (fba) = [x 0j n 



w 



= X 



(fba) 



N we find 

m ■ Xq • U N ] + [N, n M ■ AT* 

^y(fba)_ 



(37) 

ITv] 

(38) 



/,■ 



}) 



ten ' fc 

0{Tr(i? nn ) 2 } 
(34) 



/or 2-D Standard Tensor-ESPRIT^ The vector rl anaf f/ze 
matrices W ma t and Wtcn are given by 



J 2 /e JMfc 



J? 5 



1 ■ vj) ® (t/ n ■ u*) 



= (s 



S ] T 



Ti 



Im 2 ® * 



U 



>]* vH 



E4 n] c4 nlH 

-Tmi 



(35) 

A/ 2 x (Mj ■ 

(36) 



However, the latter relation shows that we are interested in the 
perturbation of the subspace of a matrix x[ ) fba ' 1 superimposed 
by an additive perturbation 

^(fba) 

which is small. Since 
the explicit perturbation expansion we have used up to this 
point requires no additional assumptions, the surprisingly 
simple answer is that we do not need to change anything 
but we can apply the previous results directly. All we need 
to do is to replace all exact (noise-free) subspaces of Xq 
by the corresponding subspaces of X Q fba K From OH), we 
immediately obtain the following explicit first-order expansion 
which is valid for R-D Unitary ESPRIT 



JV) 



A/4" 



=i m {„r )T ■ ( 



J-2 'IK 



(39) 



i Im 2 , T2 = I 



Ml 



IMi 



C2,M 2 



• AU[ iha 



(fba 

9fc 



} + 0{A 2 } 



where Af7^ fba) is given by 



and t r _ rn is the m-th column of T r . Finally, K p>q is the 
commutation matrix from (O. 

Proof: cf. Appendix Icl 

Note that for 1-D Standard ESPRIT, this MSE expression 
agrees with the one shown in lllin . However, |[T8"1 does not 
directly generalize to the tensor case. This is the advantage 
of the MSE expressions (l33l l and ( f34l > where we only need to 
replace W ma t by Wtcn to account for the enhanced signal 
subspace estimate. Furthermore, note that the special case of 
circularly symmetric white noise corresponds to R nn = a\ ■ 
Imn and C nn = Omnxmn- 

D. Incorporation of Forward-Backward-Averaging 

So far we have shown the explicit expansion and the MSE 
expressions for R-D Standard ESPRIT and R-D Standard 
Tensor-ESPRIT. In order to extend these results to Unitary- 
ESPRIT-type algorithms we need to incorporate the mandatory 
preprocessing for Unitary ESPRIT which is given by Forward- 
Backward-Averaging. The second step in Unitary ESPRIT is 
the transformation on the real-valued domain. However, it can 

4 The reason that this result is specific for R = 2 is not )34t (which applies 
to arbitrary R) but (36) which we develop only for R = 2 here. 



A[/ (fba 



U 



(fba 



u 



(fba) 



• N 



(fba) 



V 



(fba) 



,(fba)" 



(40) 



and C/< fba) , U^K V^K £( fba) correspond to the signal 
subspace, the noise subspace, the row space, and the singular 



values of X 



(fba) 




respectively. Likewise, qf^^ 



and pi ' 



represent the corresponding versions of q k and p k if U B is 
replaced by £7( fba ) in the shift invariance equations. 

With the same reasoning, an explicit expansion for R-D 
Unitary Tensor-ESPRIT is obtained by consistently replacing 
in ( |32| ). i.e., 

[/(fba)" " 



X by Af a) 



/ (fba) T / ? ( 



f(r) 



f(r) 
J n 



AW [sl(fba) 



(R+l) 



(fba) 



2 /A 
0{A 2 } 



M 

k 



f(r) 
J 1 



(41) 



Similarly, Theorem [5] can be applied to compute the MSE 
since we only assumed the noise to be zero mean and possess 
finite second order moments, which is still true after forward- 
backward averaging. The following theorem summarizes the 
results for R-D Unitary ESPRIT and 2-D Unitary Tensor- 
ESPRIT: 

Theorem 4. For the case where N or Af contain zero 
mean random variables with finite second-order moments 
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n ■ n H J and 
■ n } for 



described by the covariance matrix R nn = E j 
the complementary covariance matrix C nn = E { 

n = vcc{AT} = vccj [A/"] ^ 3) j, the MSE for R-D Unitary 
ESPRIT and 2-D Unitary Tensor-ESPRIT are given by (133 b 



1 (r) (r) 

and i34\ if we replace r k by r k 



W ton by W[y, 
<?£■>. Here, r« 



W m ,t by Wgg 



and R 

<fb. 



nn as well as C nn 

(fba) 

n n /~1 1/1/ 

(fba) 



by H^ a) and 

and W^ten** 1 flre computed 
as in (I33l l an<i ( 134b fry consistently replacing all quantities 
by their forward-backward-averaged equivalents. Moreover, 
R n ^ a ' and C^ a ' represent the covariance and the pseudo- 
covariance matrix of the forward-backward averaged noise, 
which are given by 



R 



fba 
nn 



•rffba 
'nn 



Cnn 

H-MN ■ Rn 



Cnn ■ H-MN 
n^/AT • Rtn ' HjVfJV 



Rnn ' II 



n 



MN 



MN 

n 



MN 



Proof: cf. Appendix iDl 

Note that in the special case where the noise is circularly 
symmetric and white we have R^^ = c n ' ^2MN and 

It is important to note that results on Unitary ESPRIT in this 
section relate to the LS solution only. If TLS is used instead, 
the equivalence of Standard ESPRIT with Forward-Backward- 
Averaging and Unitary ESPRIT is shown in (9). 

E. Extension to other ESPRIT-type algorithms 

In a similar manner as in the previous section, other 
ESPRIT-type algorithms can be analyzed. For instance, the 
NC Standard ESPRIT and NC Unitary ESPRIT algorithm for 
strict-sense non-circular sources are based on a different kind 
of preprocessing where instead of augmenting the columns 
we augment the rows of the measurement matrix. Yet, the 
explicit first order perturbation expansion still applies since the 
result can be written as a noise-free (augmented) measurement 
matrix superimposed by a small (augmented) perturbation 
matrix. Consequently, for the explicit expansion we only need 
to consistently replace the quantities originating from the SVD 
of Xq by the corresponding quantities from the appropriately 
preprocessed measurement matrix Xq 10 ' . Likewise, the MSE 
expressions are directly applicable since we only require 
the noise to be zero mean and possess finite second-order 
moments. 

Another possible extension is to incorporate spatial smooth- 
ing. If sources are mutually coherent, preprocessing must be 
applied to the data to decorrelate the sources prior to any 
subspace-based parameter estimation scheme. Via Forward- 
Backward-Averaging, two sources can be decorrelated. How- 
ever, if more than two sources are coherent (or if FBA cannot 
be applied), additional preprocessing is needed. For spatial 
smoothing we divide the array into a number of identical 
displaced subarrays and average the spatial covariance matrix 
over these subarrays. Since the number of subarrays we 
choose is a design parameter influencing the performance, 
investigating its effect by virtue of an analytical performance 
assessment would be desirable. Note that the spatial averaging 



introduces a correlation into the noise. Therefore, the presented 
framework is particularly attractive since for the explicit ex- 
pansion, no assumptions about the noise statistics are needed. 
A further extension is the performance assessment of tensor- 
based schemes for spatial smoothing. We have introduced a 
tensor-based formulation of spatial smoothing for R-D signals 
in ifTTI . Moreover, a tensor-based spatial smoothing technique 
for 1-D damped and undamped harmonic retrieval with a 
single snapshot is shown in 0381 . The extension to multiple 
snapshots is introduced in |39l and an R-D extension is shown 
in ll37l . A major advantage of [39], 11371 is that the performance 
of the ESPRIT-type parameter estimates is almost independent 
of the choice of the subarray size. This could be verified 
by analytical results if the performance analysis is extended 
accordingly. 

F. Incorporation of Structured Least Squares (SLS) 

So far, all performance results are based on ESPRIT using 
LS, i.e., the overdetermined shift invariance equations are 
solved using LS only. However, the LS solution to the shift 
invariance equation is in general suboptimal as errors on both 
sides of the equations need to be taken into account. Even 
more so, since for overlapping subarrays, the shift invariance 
equation has a specific structure resulting in common error 
terms on both sides of the equations, this structure should 
be taken into account when solving them. This has led to the 
development of the SLS algorithm |8 |. Since it has been shown 
that the resulting ESPRIT algorithm using SLS outperforms 
ESPRIT using LS and TLS for overlapping subarrays ||8), it is 
desirable to extend our performance analysis results to SLS- 
based ESPRIT as well. 

Due to the fact that our analysis is asymptotic in the 
SNR we can make the following simplifying assumptions for 
SLS. Firstly, we consider only a single iteration, as proposed 
in 0. This is optimal for high SNRs, since the underlying 
cost function is quadratic but actually asymptotically linear 
(the quadratic term vanishes against the linear terms for high 
SNRs). Secondly, we do not consider the optional regulariza- 
tion term in SLS (i.e., we set the corresponding regularization 
parameter a to infinity) as regularization is typically not 
needed for high SNRs. 

Under these conditions we can show the following theorem: 

Theorem 5. A first order expansion of the estimation error of 
1-D Standard ESPRIT using SLS is given by 

A^.sls =Im{r£ SLS • vec{AL/ s }} + O {A 2 } (42) 
= Im K.sls • Wmat • vcc{iV}} + O {A 2 } (43) 
where W m &t is defined in (135b and rJ SLS is given by 



fe.SLS 



pI ■ (Ji ■ u s ) + ■ 
' T (Ji ■ u s ) H 

Pk ' 



Jl 



■ Fi 



w 



R,U 



w 



R,U 



Jl) + Id® (Jl - U s (Ji ■ U a ) + ■ J; 



J 1 ■ U s (J 1 ■ U s y • Ji - (I d ® J a ) 
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SLS 



I d <g> (Ji • U s ) , (* T <g> Ji) - (J d ® J 2 ) 



for k = 1, 2, . . . , d. 77ie MSEfor zero mean noise samples can 
then be computed via 



R W 



E{(A^. SLS ) 2 } = i(rf. SLS .W: 

{rJ SLS • W mat • C nn ■ WZ* ■ r fc ,sLs} ) + O {Tr (I? nn )' 

(44) 

Proof: Equation (T43b is shown in Appendix [E] Since the 
explicit expansion of (|43T > has the same form as the explicit 
expansion in d3Qt . the MSE expression (l44t is shown analo- 
gously to Theorem [3] as presented in Appendix [C] 



G. Special case: Single source 

So far we have found closed-form expressions for the first- 
order approximate MSE for different kinds of ESPRIT-type 
algorithms. As they are deterministic, they can be plotted for 
varying system parameters without performing Monte-Carlo 
simulations and one can learn from these plots under which 
conditions the performance changes how much. 

However, it would be desirable to find expressions that are 
even more insightful. The biggest disadvantage of the MSE 
expressions in their current form is that they are formulated in 
terms of the subspaces of the noise-free observation matrix 
and not in terms of the actual parameters with a physical 
significance, such as, the number of sensors or the positions 
of the sources. 

Finding such a formulation in the general case seems to be 
impossible given the complicated algebraic nature in which 
the MSE expressions depend on the physical parameters. 
However, it becomes much easier if some special cases are 
considered. Therefore we present one example of such a 
special case in this section, namely, the case of a single 
source captured by a uniform linear array (ULA) and a 
uniform rectangular array (URA) and circularly symmetric 
white noise. Although this is a very trivial case, it serves 
as an example which types of insights such an analytical 
performance assessment can provide. For the 1-D case we 
have the following theorem: 

Theorem 6. For the case of an M -element ULA (1-D) and 
a single source (d = 1) we can show that the mean square 
estimation error of the spatial frequency for Standard ESPRIT 
and for Unitary ESPRIT is given by 



>{(Am) 2 } = - 



p (M-lf 



Oil 



(45) 



Moreover, the deterministic Cramer-Rao Bound can be sim 
plified into 

1 6 



CRB 



p M ■ (M 2 - 1) 

Consequently, the asymptotic efficiency is given by 

CRB 6(M - 1) 

71 ~ p"^o E{(Ap) 2 } ~ M(M + 1) ' 



(46) 



(47) 



Here, p represents the effective SNR given by p — where 

.Ft is the empirical transmit power given by Pt = \\S\\p /N 
if S is the matrix of source symbols. 

Proof: cf. Appendix |F| 

Note that [28] provide an MSE expression for ESPRIT for 
-.the case of a single source which scales with 1/M 2 and is 
jderived under the assumption of high "array SNR" P ■ M/a 2 , 
i.e., it is asymptotic also in M. The result presented here is 
accurate for small values of M as well and only asymptotic 
in the effective SNR N ■ -Pt/ct 2 . Also note that analytical 
expression for the stochastic Cramer-Rao Bound for one and 
two sources are available in ll33l . 

We can simplify the MSE expression for ESPRIT using 
SLS shown in Section IV-FI in a similar manner, as shown in 
the following theorem: 

Theorem 7. The mean square estimation error of the spatial 
frequency for Standard ESPRIT using SLS on an M -element 
ULA is given by 

E{(A^) 2 } 6 A/4 ~~ + 24M2 ~ 22M + 23 



P 



M(M 2 + 11) 2 (M -If 



O 



(48) 



Consequently, the asymptotic efficiency is given by 
(M 2 + 11) 2 (A/- 1) 



(M + 1)(M 4 - 2M 3 + 24Af 2 - 22M + 23) 
M 5 - M 4 + 22M 3 - 22M 2 + 121M - 121 



(49) 



M 5 - M 4 + 22M 3 + 2A/ 2 + 717 + 23 
Proof: cf. Appendix iGl 

Finally, for the 2-D case, we have the following result: 

Theorem 8. For a uniform rectangular array of Mi x Mi 
sensors (2-D) and a single source, the mean square estimation 
error of the spatial frequency for 2-D Standard ESPRIT, 2- 
D Standard Tensor-ESPRIT, 2-D Unitary ESPRIT, and 2-D 
Unitary Tensor-ESPRIT can be simplified into 



;{(A,i( 1 )) 2 + (A / i( 2 )) 2 } 



(50) 



1 

P 



1 



1 



(Mi - l) 2 ■ M 2 Mi ■ (M 2 - If 



Finally, the deterministic Cramer-Rao Bound for a URA can 
be written as 



CRB = Tr (C) 



1 

9 



6 



6 



M • (Mf - 1) M • (M| - i; 



(51) 

Proof: in Appendix [H] we derive MSE expressions for R- 
D Standard ESPRIT, i?-D Unitary ESPRIT, and the Cramer- 
Rao Bound for the more general R-D case. From these, this 
theorem follows by setting R = 2. Moreover, we show the 
identity of R-D Standard Tensor-ESPRIT and R-D Unitary 
Tensor-ESPRIT with R-D Standard ESPRIT for R = 2. 

These MSE expressions provide some interesting insights. 
Firstly, they show that for a single source there is neither 
an improvement in terms of the estimation accuracy from 
applying Forward-Backward- Averaging nor from the HOSVD- 
based subspace estimate. This is surprising at first sight since 
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the HOSVD-based subspace estimate itself is more accurate 
also for a single source. 

Moreover, they show that the asymptotic efficiency can be 
explicitly computed and that it is only a function of the array 
geometry, i.e., the number of sensors in the array. Unfor- 
tunately, the outcome of this analysis is that ESPRIT-type 
algorithms using LS are asymptotically efficient for M = 2, 3 
in the 1-D case and Mi <E [2, 3], M 2 <G [2,3] in the 2-D case. 
However, they become less and less efficient when the number 
of sensors grows, in fact, for M ->oowe even have 77 — > 0. A 
possible explanation for this phenomenon could be that an M- 
sensor ULA offers not only the one shift invariance used in LS 
(the first and last M — 1 sensors) but multiple invariances l36l . 
which are not fully exploited by LS. However, for ESPRIT 
based on SLS, the asymptotic efficiency is in general higher, 
in fact, we have r) — >• 1 for M = 2, M = 3 and M 00 for 
a single source. Moreover, even for limited M, r] is never far 
away from 1. As we show in the simulations below, we have 
77 = 1 for M = 2,3 and the smallest value of 77 is obtained 
for M = 5 where 77 = 36/37 0.973 for d = 1. 



B 

CO 



CD 

o 10 

nS 
Q_ 
u) 

</> 

o 

LU 
CO 



10' 



25 



30 



s ^"l s V_ 


-V" Matrix: Empirical 
-©-Matrix: Analytical [n] 
-^-Matrix: Analytical [n] and [s] 
-A- Tensor: Empirical 
-B- Tensor: Analytical [n] 

Tensor: Analytical [n] and [s] 


^^^s 





35 40 
SNR [dB] 



45 



50 



0.7, fi> 



and 



Figure 1. Subspace estimation accuracy using only vs. using T 
T^ s >. Scenario: d = 3 correlated sources (p = 0.97) at //" 

0.9, /4 1 ' = l.l,Mi 8) = -0.1,At2 2) = -0-3, ^ = -0.5, an 8 x 8 URA 
and N = 20 snapshots. 



VI. Simulation results 

In this section we show numerical results to demonstrate the 
asymptotic behavior of the analytical performance assessment 
presented in this chapter. We first investigate the subspace 
estimation accuracy in order to verify d26| i. Note that the 
analytical results for the subspace estimates are explicit ex- 
pansions in terms of the perturbation (i.e., the additive noise). 
Therefore, we repeat the experiment with a number of ran- 
domly generated realizations of the noise and perform Monte- 
Carlo averaging over the analytical expansions. These "semi- 
analytical" results are then compared with purely empirical 
results where we estimate the subspace via an SVD or a 
HOSVD and compute the estimation error compared to the 
true signal subspace. 

The subsequent numerical results demonstrate the perfor- 
mance of ESPRIT-type parameter estimation schemes. Here, 
we compute the mean square estimation error in three different 
ways. Firstly, analytically, via the MSE expressions provided 
in Theorem [3] Theorem |H and Theorem [5] (eqn. d44b ~). re- 
spectively. Secondly, semi-analytically, by performing Monte- 
Carlo averaging over the explicit first-order expansions pro- 
vided in equation d32b . d4Tb . and d43t . respectively. Thirdly, 
empirically, by estimating the spatial frequencies via the 
corresponding ESPRIT-type algorithms and comparing the 
estimates to the true spatial frequencies. 

For all the simulations we assume a known number of planar 
wavefronts impinging on an antenna array of M istrotropic 
sensor elements. We assume uniform A/2 spacing in all 
dimensions, i.e., an M-element uniform linear array (ULA) in 
the 1-D case and an Mi x M2 uniform rectangular array (URA) 
in the 2-D case. The sources emit narrow-band waveforms 
Si(t) modeled as complex Gaussian distributed symbols Si(i) 
and we observe N subsequent snapshots t = 1, 2, . . . , N. All 
sources are assumed to have unit power, i.e., E {|si(£)| 2 } = 1. 
In the case where source correlation is investigated we gener- 
ate the symbols Si(t) such that E {s,(i) ■ Sj(t)*} = p-eP^ for 



i j = 1,2, ... ,d, where p is the correlation coefficient be- 
tween each pair of sources and ipij is a uniformly distributed 
correlation phase. The additive noise is generated according 
to a circularly symmetric complex Gaussian distribution with 
zero mean and variance and noise samples are assumed to 
be mutually independent. Therefore, the Signal to Noise Ratio 
(SNR) is defined as 1/al 

A. Subspace estimation accuracy 

We evaluate the subspace estimation accuarcy by computing 
the Frobenius norm of the subspace estimation error, i.e., 



Ai7 s II u in the matrix case and 



AW 



(R+l) 



the 



tensor case. 

In order to find the estimation error empirically, we obtain a 
subspace estimate (7 S via an SVD of the noisy observation and 
then compare it to the true subspace U s column by column. 
The estimation error of the n-th column is computed via 



Au r , 



1,2, 



(52) 



to account for the inherent phase ambiguity in each column 
of the SVD, cf. f!9l . 

For the analytical estimation error we calculate At7 s via 
the first-order expansion AU S ~ U n ■ r' n ' provided in (fT9l 
and the expansion AU S U n ■ T [n] + U s ■ T [s] provided 
in ( 1201 , respectively. Note that the latter is more accurate since 
it additionally considers the perturbation of the individual 
singular vectors, i.e., the particular choice of the basis for the 
signal subspace. However, this contribution is irrelevant for 
the performance of ESPRIT-type algorithms. 

In Figure [TJ we have d = 3 sources positioned at p^ = 

1.9, /4 X) = i-i, Ati 2) = -o.i,/4 2) = -°- 3 > /4 2) = 

0.5 and mutually correlated with a correlation coefficient of 



0.7, p<£ ] 
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Figure 2. Subspace estimation accuracy using ]S n ' only vs. using Y 

(1) - -1.5, & 
0.2, /4 2) = 0.7, fj,f 



rM. Scenario: d 
— 1.5, an 8 X 



1.0, 



4 uncorrelated sources at ^ 

-0.3, pf ) = 1.3, } 
URA, and N = 5 snapshots. 



and 

0.5, 



Figure 3. Performance of 2-D SE, STE, UE, UTE for d 
sources (p = 0.9999) located at p^ 1 ' = 1, 
/4 2) = 1, a 5 X 6 URA, and AT = 20 snapshots 



2 highly correlated 
0.5, and 



0.5, p< 2) 



p = 0.97. Moreover, the array size is increased to an 8 x 8 
URA. For the simulation result shown in Figure [2] we consider 



d = 4 uncorrelated sources located at p^ = —1.5, fin 1 ' = 0.5, 
pi ] = 1.0, p^ = -0-3, a4 2) = !-3. M2 2) = -0-2. H ] = 0.7, 
/if ) = -1.5 and iV 



-0.3, /xi' 

5 snapshots. 

Both simulations show that the empirical estimation errors 
agree with the analytical results as the SNR tends to infinity. 
Therefore, the improvement obtained by the HOSVD-based 
subspace estimate can reliably be predicted via the analytical 
expressions. In general, it is particularly pronounced for corre- 
lated sources and for a small number of snapshots. Moreover, 
while for three correlated sources, the impact of the additional 
term U s ■ is negligibly small, it is clearly visible for the 
four uncorrelated sources shown in Figure [2] 

B. R-D Tensor-ESPRIT 

The following set of simulation results demonstrates the 
performance of R-D matrix-based and tensor-based ESPRIT. 
As explained in the beginning of this section, for the analytical 
results we use Q and (O for R-D Standard ESPRIT and R- 
D Standard Tensor-ESPRIT, respectively, and their extensions 
for Unitary ESPRIT and R-D Unitary Tensor-ESPRIT as 
discussed in Theorem |U Likewise, the semi-analytical results 
are obtained by Monte-Carlo averaging of the explicit first- 
order expansion provided in d32b and fiTJ ). respectively. 

For Figure [3] we employ a 5 x 6 URA and collect N = 



20 snapshots from two sources located at p^' = 1, p\ 
—0.5, p± = —0.5, and p£ — !■ The sources are highly 
correlated with a correlation of p = 0.9999. On the other hand, 
for Figure [4] we increase the number of sources to d = 3 and 
the correlation coefficient to p = 0.97. Moreover, the spatial 



frequencies of the sources are given by p^ - 

0.9, $ = LI, /4 2) = -0.1,a4 2) = -0.3,a4 2) 
we use an 8 x 8 URA. 



0.7, pf 



-0.5 and 
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Figure 4. Performance of 2-D SE, STE, UE, UTE for d = 3 correlated 



sources (p 

i.i,/4 2) = 

snapshots. 



0.97) positioned at p 



(i) 



-o.i,/4 2) 



-0.3, /4 2) 



0.7, 



0.9, 



-0.5, an 8 X 8 URA, and N = 20 



To enhance the legibility, we show the semi-analytical 
estimation errors only in Figure [3] since they always agree with 
the analytical results, as expected. Moreover, the empirical 
estimation errors agree with the analytical ones for high 
SNRs. This is also expected as the performance analysis 
framework presented here is asymptotically accurate for high 
effective SNRs. We conclude that the improvement in terms 
of estimation accuracy for Tensor-ESPRIT-type parameter 
estimation schemes can reliably be predicted via the analytical 
expressions we have derived. 
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Figure 5. Performance of LS-ESPRIT vs. SLS-ESPRIT for 4 sources at 
HI = 1.0, (U2 = 0.7, i^s = -0.6, (M = -0.3, an M = 8 ULA, N = 3 
shapshots. 



Figure 7. Performance of LS-ESPRIT and SLS-ESPRIT for a single source 
vs. the number of sensors M (M-ULA) at an effective SNR of 25 dB (P T = 
1, 0-2 = 0.032, N = 10). 
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Figure 6. Performance of LS-ESPRIT vs. SLS-ESPRIT for d = 3 correlated 
sources (p = 0.99) at fii = 1, fi2 = 0, [13 = -1, a M = 12 ULA and 
N = 10 shapshots. 



C. Structured Least Squares 

The next set of simulation results illustrates the analytical 
expressions for ESPRIT using SLS. The semi-analytical MSE 
is obtained by Monte-Carlo averaging over the explicit ex- 
pansion provided in d43l > and the analytical MSE is computed 
via ( 144-b . For the empirical estimation errors we perform a 
single iteration of the Structured Least Squares algorithm and 
do not use regularization (i.e., the regularization parameter a 
is set to oo). 

The first simulation result is shown in Figure [5] Here we 
consider N = 3 snapshots from d = 4 uncorrected sources 
captured by an M — 8 element uniform linear array. The 
sources' spatial frequencies are given by \i\ = 1.0, fj,2 = 
0.7, 113 = —0.6, /14 = —0.3. Note that since N < d, we cannot 



apply Standard ESPRIT, therefore, only Unitary ESPRIT is 
used. On the other hand, in the second scenario we consider 
N = 10 snapshots from d = 3 sources that are mutually 
correlated with a correlation coefficient of p = 0.99. The 
sources are located at \x\ = 1, ^2 = 0, (13 = —1 and a 
M = 12 element ULA is used. The corresponding estimation 
errors are shown in Figure [6] 

As before, the empirical results agree with the analytical 
results for high SNRs. Moreover, the improvement in MSE 
obtained via SLS is particularly pronounced for the correlated 
sources. However, even the very slight improvement which is 
present for four uncorrelated sources can reliably be predicted 
via the analytical MSE expressions we have derived. 

D. Asymptotic efficiency for a single source 

The final set of simulation results demonstrates the special 
case of a single source, in which case the MSE expressions can 
be simplified to very compact closed-form expressions which 
only depend on the physical parameters, i.e., the array size 
and the SNR. 

Figures [7] and [8] show the MSE vs. the number of sensors 
M for a (1-D) Uniform Linear Array and vs. M\ for a (2-D) 
Mi x Mi Uniform Rectangular Array, respectively. For both 
scenarios, the spatial frequencies of the single source were 
drawn randomly (note that they have no impact on the MSE). 
The effective SNR was set to 25 dB for Figure Q (P? = 
1,0-2 = 0.032, N = 10) and to 46 dB for Figure E (P T = 
1, al = 10~ 4 , N = 4), respectively. 

For both plots we observe that LS-ESPRIT is asymptotically 
efficient for M = 2 and M = 3 (which, in the 2-D case means, 
a 3 x 3 URA) and then becomes increasingly inefficient as the 
array size grows. Moreover, for the 1-D case we see that SLS- 
based ESPRIT is in fact very close to the Cramer-Rao Bound, 
which may, at first sight, lead to believe that the asymptotic 
efficiency is in fact 1 for all M. However, as we have shown 
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Figure 9. Asymptotic efficiency of LS-ESPRIT vs. SLS-ESPRIT vs. M. Same scenario as in Figure [7] The left-hand side shows a zoom. 
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Figure 8. Performance of LS-ESPRIT for a single source vs. Mi using an 
Ml X Mi URA at an effective SNR of 46 dB (Pp = 1, <rl = 1CT 4 , N = 4). 

it is in fact slightly lower than one. Therefore, we provide two 
additional figures where we depict the "asymptotic efficiency", 
i.e., we divide the CRB by the corresponding value of the 
MSE. The resulting efficiency plot is shown in Figure [9] This 
plot shows more clearly that LS-ESPRIT becomes increasingly 
inefficient for M > 3, whereas SLS-ESPRIT approaches rj = 
1 for large M, The worst efficiency is found for M = 5 where 
we have rj = 36/37 « 0.973. 

VII. Conclusions 

In this paper we have discussed a framework for analytical 
performance assessment of subspace-based parameter estima- 
tion schemes. It is based on earlier results on an explicit 
first-order expansion of the SVD and its application to 1- 
D versions of subspace-based parameter estimation schemes, 
e.g., ESPRIT. We have extended this framework in a number of 
ways. Firstly, we have derived an explicit first-order expansion 
of the HOSVD-based subspace estimate which is the basis 
for Tensor-ESPRIT-type algorithms. Secondly, we have shown 
that the first-order expansion for 1-D Standard ESPRIT can 



be extended to other ESPRIT-type algorithms, e.g., R-D Stan- 
dard ESPRIT, R-D Unitary ESPRIT, R-D Standard Tensor- 
ESPRIT, or R-D Unitary Tensor-ESPRITThirdly, we have 
derived a corresponding first-order expansion for Structured 
Least Squared (SLS)-based ESPRIT-type algorithms. 

All these expansions have in common that they are explicit, 
i.e., no assumption about the statistics of either desired signal 
or additive perturbation need to be made. We only require the 
perturbation to be small compared to the desired signal. We 
also do not need the number of snapshots to be large, i.e., 
they even apply to the single snapshot case (N = 1). Our 
fourth contribution is to show that the mean square error can 
readily be computed in closed-form and that it depends only on 
the second-order moments of the noise. Consequently, for the 
MSE expressions we only need the noise to be zero mean and 
its second order moments to be finite. Neither Gaussianity nor 
circular symmetry is required. This is a particularly attractive 
feature of our approach with respect to different types of 
preprocessing which alters the noise statistics, e.g., spatial 
smoothing (which yields spatially correlated noise) or forward- 
backward averaging (which annihilates the circular symmetry 
of the noise). Since we do not require spatial whiteness or 
circular symmetry, our MSE expressions are directly appli- 
cable. The resulting MSE expressions are asymptotic in the 
effective SNR, i.e., they become accurate as either the noise 
variance goes to zero or the number of observations goes to 
infinity. 

As a final contribution we have investigated the special 
case of a single source, circularly symmetric white noise, and 
uniform linear (1-D) or uniform rectangular (2-D) arrays. In 
this case we have been able to show analytically, that R-D 
Standard ESPRIT, R-D Unitary ESPRIT, and (for R = 2) R- 
D Standard Tensor-ESPRIT as well as R-D Unitary Tensor- 
ESPRIT yield the same MSE, which only depends on the 
effective SNR and the number of antenna elements. We have 
also shown that 1-D Standard ESPRIT using SLS has a lower 
MSE which is also expressed explicitly as a function of the 
effective SNR and the number of antenna elements. 
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Appendix A 
Proof of TheoremQ] 

As shown in ( [Tol l, the estimated signal subspace tensor can 
be computed via 



U 



.a.— 1 



S x x U 1 ... x R U R x r+1 S s 



(53) 



Here, <S represents the truncated version core tensor S from 
the HOSVD of X. In order to eliminate <s' ' we require the 
following Lemma: 

Lemma 1. The truncated core tensor <s' ' can be computed 
from X directly via 



X X! U[ 



TT l ' 

Xr+i Ur+i- 



(54) 



Proof: To show d54l i we insert the HOSVD of X given 
by X = S Xi Ui . . . X R +i U R +i. Using we obtain 



S = Sx l [U[' -Ui 



Xfl+i [U R+1 •Ur+x). (55) 



However, since the matrices of r-mode singular vectors U r 

- [s] H 

are unitary, they satisfy U r ■ U r = [l Pr , PrX (M r -p r )\ ■ 
Therefore, <s' ' computed via ( l54i i contains the first p r ele- 
ments of S in the r-th mode, which shows that it is indeed 
the truncated core tensor. ■ 
Next, we use Lemma [T] to eliminate S in (1531 . We obtain 



(56) 



U =X Xl [U[' -U\' ...Xr[U-.U r 

XR+1 I ^ s ' Ur+i 

=X Xifx. 



x r T R x R+1 [ S = 1 



R+l I ' 



(57) 



where we have introduced the short hand notation T r = U r 

-m h r~[sr T 



U r . The next step is to compute the matrix 
Inserting (l57l i and using ©, we obtain 



(R+l) 



- 8 



^ = (T ! ® . . . ® f fi ) ■ • E/^ • £ s \ 

(58) 



As pointed out in Section IIII-AI the link between the 
measurement matrix X and the measurement tensor X is 
given by X = [X]T R+1) . Therefore, their SVDs (cf. (HD 
and d23l ) are linked through the following identities 

l> -t/ [sl * rV -v ln] * v -u [s] * v -r/ [nl * 

U s — V R+l, >J n — V R+1 , V s ~ U R+1 , V n — U R +i- 

Consequently we can write 



[X]Jr+1) ' ^fi+l ' — ^ ' ' S s 



Finally, inserting d59l l into 

J (ii+i) 

which is the desired result 



yields 



Ti 



..®T R )-U 



(59) 

(60) 
□ 



Corollary 1. A corollary which follows from this theorem is 
that the exact subspace U s satisfies the following identity 



U s = (Ti 



T R ) ■ U s 



(61) 



Proof: The corollary follows by considering the special 

U s . 



case where X = Xn and hence T r = T r as well as U 

T r ,,iT 



For this case we also have 



U 



(R+l) 



(R+l) 



u s , 



where the last identity resembles the fact that in the noise-free 
case, the HOSVD-based subspace estimate coincides with the 
SVD-bases subspace estimate. ■ 

Appendix B 
Proof of Theorem[2] 

We start by inserting U B = U S +AU S and T r = T r + AT r 
into d25l l. Then we obtain 

T 



u 



(R+l) 



= [(Ti + ATi) ® . . . ® (T R + AT R )} ■ {U s + AU S ) 



Us- 



[ATi ® T 2 

[Ti ® T 2 <g> 



®T R ]-U B + ... 
AT R ]-U s +0{A 2 } , 



T R ] ■ AU S 
(62) 



since all terms that contain more than one perturbation term 
can be absorbed into O {A 2 }. The first term in (l62l l represents 
the exact signal subspace (cf. Corollary [TJ, hence the remain- 



AU [H 



As 



I -R+i) 



ing terms are the first order expansion of 

the first term of this expansion already agrees with Theorem[2] 
we still need to show that for the remaining terms we have 
for r = 1,2,..., R 



Tit 
Ti 



AT r ®...T R ]-U s 
) (f/W ■ rfi • U [ f) C 



U a 



(63) 
0{A 2 }. 



As a first step, we expand the left-hand side of d63l by 
applying Corollary Q] 

[Ti ® . . . ® AT r ®...T R ]-U B 
= [Ti ® . . . <g> AT r ® . . . T fl ] ■ [Ti ® . . . ® T r ® . . . Tr] ■ l/ s 
= [(Ti • Ti) ® . . . g> (AT r • T r ) »...(T R - T R )} ■ U s 
= [Ti ® . . . <g) (AT r ■ T r ) <g> . . . r fl ] • t/ s , (64) 

where we have used the fact that the matrices T r are projection 
matrices and hence idempotent, i.e., T r ■ T r = T r . What 
remains to be shown is that AT r ■ T r = U l " ] ■ ■ U [ ° ]H + 

0{A 2 }. Since f T = U [ * ] ■ U [ f and U^ = + AU®, 
a first order expansion for AT r is obtained via 

[/[?!+ A[/M) • (u^ H + AU^ H ) 

T r + U^ ■ AU^ ]H + AU [ ° ] ■ U^ ]H + O {A 2 } 
Uf ■ AU [ f + AU^ ■ U^ + O { A 2 } , (65) 



AT r 



where in general we have AUf' = U r n ' ■ T 
O {A 2 } (cf. d24]i). Using this expansion in 

ATr ■ T r 



up ■ rf - 

we obtain 
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=Lf ] ■ (cf nl -if! + Lf ] -if) -T r 
+ (u r n] ■ if 1 + Uf ■ if 1 ) • Lf ]H • T r + O { A 2 } 

=t/|f] ■ if i H • u^ H ■ T r + t/f i ■ if H • up" ■ T r 
+ uf ■ if • Lf ]H ■ T r + u r n] • if 1 • Lf ]H ■ t t + o {a 2 } 
=Lf i ■ if i H ■ u^j^+u® ■ (if + if H ) -U^f ■ T r 

O.M,. ;> r .x M T n 



t/|. n] • if 1 • Lf ]H • T r +0{A 2 } 



-ui 



U [ * ] +0{A 2 }, 



(66) 







which is the desired result. Note that Tj. J + ly J 
follows from the fact that if' is a skew-Hermitian matrix 
(which is apparent from its definition shown in (12410. This 
completes the proof of the theorem. □ 

Appendix C 
Proof of Theorem[3] 

For R-D Standard ESPRIT, the explicit first-order expansion 
of the estimation error for the A;-th spatial frequency in the r- 
th mode in terms of the signal subspace estimation error AU S 
is given by ( f30b . This error can be expressed in terms of the 
perturbation (noise) matrix TV by inserting (fT9b . We obtain 

A^ r) =Im {r<: r)T • vec {AU S }} + O {A 2 } 

=Im {r> )T ■ W mat ■ vec {TV}} + O {A 2 } (67) 



-W 1 

' k 

w mat = 



(r) J 
Ik 



s- 1 • v T 



^2 V^i r) _ J\ 



which follows directly by applying property ([T]i to ( f30b and 
to (O. In order to expand E j(A/^) 2 j using (|67]i, we 
observe that for arbitrary complex vectors Zi,z 2 we have 
lm{zj ■ z 2 } = Im{zi} T ■ Rc{z 2 } + Re{zi} T ■ Im{z 2 } 



and hence 



Im{zJ • z 2 y 



Im{zi} T ■ Re{z 2 } ■ Rc{2 2 } T ■ Im{zi} 

+ Rc{z 1 } T -Im{z 2 } ■Im{^ 2 } T - Re {2:1} 

+ Im{2 1 } T -Rc{z 2 } ■Im{2 2 } T - Re {2:1} 

+ Re { Zl } T • Im {z 2 } • Re {z 2 } T ■ Im {zi} 

(68) 

Using ([67) in E |(A J u^' ) ) 2 | and applying Q for zj = 
rl ■ W ma t and z 2 = vec {N} = n we find 
E j(A^ r) ) 2 j =E|lm{z^} •Rc{n}-Re{n} T -Im{2i}| 
+ E|Re{z^} ■Im{n}-Im{n} T -Rc{2i}| 
+ E|lm{zJ} • Re{n} • Im{n} T • Re{zi}} 

+ E jRc{2^} ■ Im{n} ■ Rc{n} T • Im{zi}} 

(69) 



Since the only random quantity in (1691 is the vector of noise 
samples n, we can move Z\ out of the expectation operator. 
We are then left with the covariance matrices of the real part 
and the imaginary part of the noise, respectively, as well as 
with the cross-covariance matrix between the real and the 
imaginary part. To proceed we require the following lemma: 

Lemma 2. Let n be a zero mean random vector with co- 
variance matrix R nn = E{n-n H } and pseudo-covariance 
matrix C nn — E jn • n T }. Then, the covariance matrices of 
the real part of n, the imaginary part of n and the cross- 
covariance between the real and the imaginary part of n are 



given by 










= E 


^Re{n} 


Re{n} T | = 


^Rc{R Dn + 


C nn } 


= E 


Im {n} 


Im{n} T | = 


-Re {iinn — 


Cnn} 


= e 


^Re{n} 


Im{n} T | = 


--Im{-Rnn 


— C n n} 


R^ = E 


j^Im {n} 


Rc{n} T | = 


-Im{-R nn + 


C nn } ■ 



Proof: To prove this Lemma we expand R nn and C„ 
by inserting n = Re {n} + jlm {n}. We then obtain 



Rnn = R^ + + J (l& R^) (70) 

Cnn = R^ R) - ^ + J + R^) ■ (71) 

Since R D ^ R) , R^^, R ( nf \ and Rnn ] are real-valued, the 
solution of ( l70l i and (T7TT > is straightforward. ■ 
Using Lemma |2] in (f69t we obtain for the mean square error 

E{(A^ r) ) 2 } = l(hn{zi} ■Rc{ J Rnn + C n „}-Im{z 1 } 
+ Re {zj} ■ Re {R nn - C n „} • Re { Zl } 

+ Im {zj} ■ lm{-R nn + Cnn} ■ Rcj^l} 

+ Re {zl} ■ Im{R n „ + C n „} ■ Im{zi}) 

(72) 



Finally, ( |72l can be expressed in more compact form as 



E 



{(A M M) 2 } = l(zf- R^ ■z 1 -Rc{zJ- Cl n ■ Zl }) 

(73) 



for Z\ = W^„+ • Tu, which is the desired result 



□ 



mat k 

Note that Gaussianity is not needed for these properties to 
hold. Consequently, the MSE expressions are still valid if the 
noise is not Gaussian. We only need it to be zero mean. Also, 
note that for the special case of circularly symmetric white 
noise we have R nn = o\ ■ Imn an d C nn = OmnxMN and 
hence the MSE simplifies into 



■ W: r) ) 2 } = f 



n 2 
l*ill 2 



2 



Wl 



,(r) 



(74) 



The procedure for 2-D Standard Tensor-ESPRIT is in fact 
quite similar. The first step is to express the estimation 
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error in fj,^ in terms of the perturbation n = vec {TV} 
vcc {W}J 3} } (cf. CI)- This expression takes the form 



An 



(r) _ 



Im 

Im 

„i-,T 



vec 



AiT 



(R+l) 



Uy ■ up 

T 2 - 



Uj ® Im 



f M 
| r fc 



-0{A 2 } 

W ton -vec{TV}}+0{A 2 } (75) 



vec ■ 



{ [u [ * ] ■ U [ ? 



m (1) 



t/J (8 Im 



AW 



depends linearly on vcc {TV}. Due to the 

J (R+l) 

fact that d75l l has the same form as (IFTT i, the second step to 
expand the MSE expressions follows the same lines as for R- 
D Standard ESPRIT, which immediately shows that the MSE 
becomes 

:{(A M W) 2 } = !(*?. h* . Zl Rc {,t . c x . zi }) 

(76) 



T 2 - 



®T 2 } 
} (80) 



where the matrix T 2 is constructed from the columns of T 2 
given by t 2 , m for m = 1,2,..., M 2 in the following manner 



E 



■Mi 



Imi ® *2,1 



- Mi 



(81) 



for zi = VF^, n • r£ . Therefore, the final step is finding an 
explicit expression for Wten which satisfies 

T 



AW lsl 



(3) 



vec{TV} + 0{A 2 }. (77) 



The final step is to rearrange the elements of vcc i [A/] ^ j so 
that they appear in the same order as in vcc {TV}. However, 
since TV = [A/] (3) , this can easily be achieved in the following 



manner 



Recall from Theorem that for R = 2, the HOSVD- 

r . [ s ] "I 

based signal subspace estimation error AU can be 

expanded into 



= K M2X{Ml . N) ■ vec {AT} 



(82) 



(R+l) 



AZT 



(:s) 



(Ti (» T 2 ) • At/ S 



AC/ s| • E7 



*7 



where ^M 2 x(Mi-iV) is the commutation matrix (cf. equa- 
tion (f2ll). This completes the derivation of the second term 
of Wt en . The third term is obtained in a similar manner. In 



T 1 



au 1 ? ■ m 



U s + 0{A 2 }, 



(78) 



this case, no permutation is needed, since vec 
vcc {[A/iy = vcc {TV} . 



□ 



where AU S , AU\ , and AU 2 are given by 
AU S = U a ■ U* • TV • V s • 



TV-C/£ 



AC/Jf 1 = t/| n] ■ C/[, n]H ■ [M] (r yV 



J :>, 

si 



and 

for r 



1,2. 
(79) 



Appendix D 
Proof of Theorem^] 

As pointed out in Section IV-DI the inclusion of Forward- 
Backward- Averaging leads to a very similar model, where all 
quantities originating from the noise-free observation Xq (or 
Xq) are replaced by the corresponding quantities for Xq 



(fba) 



The first term in (1781 is easily vectorized by applying prop 
erty <JTJ which yields the first term of Wten as 

vec{(T 1 ®T 2 )-AC/ s } 

= vec{(Ti®T 2 )- v' 11 " 



vk 



TV • (7 



= cn s| -S 



= s 



(Ti ® r 2 ) • v^ 1 -v^ 1 
(Ti®r 2 )- v 3 n] * -y' nlT 



(or A' fba ' ) ). Since for Theorem [3] it was only assumed that 
the desired signal component is superimposed by a zero mean 
noise contribution, it is directly applicable. 

The only point we need to derive are the covariance matrix 
i and the pseudo-covariance matrix of the forward-backward 

J averaged noise n^ fba ' = vec jiV^'H, which are needed for 

• vec {TV} me MSE expressions. To this end, we can express n^ fba ' as 



vec {TV} . 



vcc 



However, for the second term in (1781 we get 



Im 



[Af] (ir vf -E^ -U® 



vcc I ( 

(t/J ® I M ) • 

{ [c/i nl • [/i n]H • • *f • SF 1 • C/I S]H ] ® T 2 } 



2 )-t/s} 



vec {TV} 
(11^® n M ) • vec {TV*} 



n 

U NM ■ n* 



(83) 



vec ■ 



by inserting (|79l for AL/j . To proceed we need to rewrite 
the vectorization of a Kronecker product. After straightforward 
calculations we obtain 

(Uj ® I M 



Equation CI83D allows us to express the covariance matrix and 
the pseudo-covariance matrix of n^^' via the covariance 
matrix and the pseudo-covariance matrix of n. We obtain 



{' 



,(fba) 



n 



(fba 



"} 

E{n ■ n H } 
IIm;v -E{n* -n H } 



E{n • n 1 
rW-E{n* 



}n 



MAT 



MAT 
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fba 



Cnn 
IIfcf.IV ■ R 



E | n (fba) . n (fba) T | 



E{ 



n ■ n 



n MN ■ e \ 

C nn 



n 



Rnn 



E{n-n H } 'IImjv 
Umn ■ E { 



This completes the proof of the theorem 



□ 



Appendix E 
Proof of Theorem[5] 

Without regularization, the cost function for 1-D Structured 
Least Squares can be expressed as J8) 



SLS 



arg mm 
*, At7 s 



J x ■ [U s + AU S • * - J 



h ■ (u s 



AU S 
(84) 



where we have used AU S only to avoid confusion with the 
AU S associated to the estimation error in U s - Note that 
the cost function is solved in an iterative manner starting 
with AU S = Onixd and with = \I>ls, where \I>ls = 
^ J i ■ C/ s ) ■ ( J 2 ■ U$ \ represents the LS solution to the shift 
invariance equation. As we compute only a single iteration 
we find one update term for U s and one for 'S'ls which we 
denote as AU s ,sls and A'S'sls (i- e -> AU s .sls represents the 
AU S which minimizes the linearized version of d84b). In other 
words, the cost function becomes 



*sls =*ls + A* SLS where 



(85) 



A* 



SLS 



arg mm 
A*,a77 s 



- J 2 - (U s + AU t 



Ji (u s + AU a ) ■(* 



'LS 



A*) 



arg mm 
A* : At/ s 



R LS + Ji • AU, 
AU s + 0{A 2 } 



LS 



Ji-U s - A* 
(86) 



where we have defined the matrix i?LS = Ji • U s ■ 'S'ls — J 2 • 
U s which contains the residual error in the shift invariance 
equation after the LS fit. Since (|86] l is linearized by skipping 
the quadratic terms in O {A 2 }, it is easily solved by an LS fit. 
To express the result in closed-form we vectorize (|86i l using 



the fact that I LAI 



|vec{A}j| 2 and obtain 



A ^SLS 

Am SjS ls 



- ~ F SLS • r LS 



(89) 



where we have introduced the vectorized quantities t*ls = 
vec{i?Ls}, Ai/'sls = vcc{A* SLS }, and Am SjS ls = 
vec {At/ s sLs}> respectively. We have also skipped the 
quadratic terms C{A 2 } in ( T88T > for the solution in d89b . 
Moreover, the matrix F S ls € C (M ~ 1)dx(<i2+M ' d) becomes 

: T 



SLS 



Ji U s 



LS 



(90) 

Therefore, our next goal is to find a first order expansion of 
A-0 SLg in d89l . This looks difficult at first sight as it involves 
an expansion of a pseudo-inverse due to -Fsls- However, this 
step simplifies significantly by realizing that .Fsls can be 
expressed as F S ls = -Fsls + AF S ls, where 

2 „ 



SLS 



I d ® (Ji • U s ) , (* T <g> Ji) - (I d <g> J 2 



(91) 



AFsls = 



(Ji • AU S ) , (A* 



Is® Ji 



where \&ls = * + A^ls- Since Fsls is n °t random (i.e.. 
only dependent on Xq but not on TV) and A'S'ls = Odxd + 
O {A}, i.e., at least linear in the perturbation, we have 



SLS 



F+ LS +0{A}. 



(92) 



A</> 



sls = arg mm _ 

vcc{A*},voc{AI7 3 } 



t-ls 



LS 



Jl 



vec 



I d ®[Ji-U B )) -vec {A*} 
- [Id g> J 2) ■ vec { At7 s } + O {A 2 } 



A ^SLS 

+ F S LS 



arg mm 

vec { A * } , vcc{ AT7 S } 

vec {A*} 
vec{AC7 s } 



C{A 2 } 



{AU S } 



(87) 



(88) 



This relation only describes the "zero-th" term of the expan- 
sion of F gLS . However, as we see below, the linear term is 
not needed for a first order expansion of Ai/> SLS . Continuing 
with (|89l , the second term of the right-hand side is given by 
t*ls for which we can write 

r LS =vec | Ji ■ U s ■ * LS - J 2 • U s } 

=vec{ J 1 ■ (U s + AU S ) ■ (* + A* LS ) 

- J 2 • (U s + AU S )} 
=vec{ Ji • AU S * + Ji U B A* LS 

-J 2 -AU s +0{A 2 }} (93) 

Moreover, as shown in [17|, A'S'ls can be expressed in terms 
of AU S via 

A* LS =( Ji • U s )+ ■ J 2 ■ AUs 

- (Ji • U s )+ • Ji • AU S ■ * + O {A 2 } . (94) 

Using this expansion in d93l we obtain 

r LS = vcc{Jj ■ AU S ■ *} 

+ vec{J 1 -U a -(J 1 -U a ) + -J 2 -AU B } 
vec{Ji-[/ s -(J 1 -[/ s )+- J!-AC/ S -*} 
vec{J 2 - AU s } + 0{A 2 } 
W R ,u-vcc{A[/ s } + 0{A 2 } where (95) 

T ® Ji) + i d ® ( Ji • Us ( Ji • u s y 

T ® (jl • U a ( Jl ■ U s ) + ■ Jl) -(I d ( 

Using ([95) and ((92) in ([89) we find that 



W RiU = ( 

and 



• J 2 

J2) 



A V'SLS 

Aw SjS ls 



- ~ F S-LS ' r LS 
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SLS 



0{A}) ■ (W R ,u ■ vec{A«7 s } + O {A 2 }) A. MSE for Standard ESPRIT 

(96) We start by simplifying the MSE expression for 1-D Stan- 
dard ESPRIT. In the case of a single source we can write 



-^sls ' ^R,u ■ vcc{A(7 s } + O {A 2 } 



Note that from d9TT > it follows that .Fsls has full row-rank 
and hence its pseudo-inverse can be expressed as -F"sls = 



SLS 



(F 



SLS 



F 



H n-1 
SLSJ 



This allows us to extract Axp SLS 
from d96l l since Ah S) slS is not explicitly needed as long as 
only one SLS iteration is performed. We obtain 



A</> 



SLS 



/ d ®(J!-[/ s ) H -(Fsls-F 



SLS J 



W RiV -vec{AU s } + 0{A 2 } 



(97) 



The final step in SLS-based ESPRIT is to replace the LS- 
based estimate * LS by * SLS = *ls + A* SLS = * + 
A'S'ls + A\l/gLs- Following the first-order expansion for 
Standard ESPRIT from HD we obtain 

Afi k = Im {pj ■ (A* LS + A* S ls) ■ Q k ) h mk + O {A 2 } 
= Im{pT.A* LS -q fc }/c^ 

+ Im{(ql ® pi) • AV SLS } /V« + O { A 2 } (98) 

Since the first term is exactly the same as for LS-based 1-D 
Standard ESPRIT, we can use the result from fl7l . Inserting 
the expansion for Ai/> SLS from ( |97| i we obtain 



A Mfc =Im<!K • (Ji-l/ 8 r- ( -^-Ji )-AU s -q k 



H( 

• Wr,u ■ vec 



ql®pl-{Ji-U s 



H s-1 



(F SLS • -F S ls) 



{A17.} }/e»" +0{A 2 } 



Finally, rearranging the first term as a T ■ B ■ c = (c T i 
vec{B} we have 



(99) 
>a T )- 



Im{- 



9fc»lPfc -(Ji-U s ) + 



J 2 



Jl 



■{AU S }} 



Im 



{( 



T „ T (" 1 ' U s 

qt ®pl ■ 



H \-l 



(F SLS • -Fsls) 



fc,SLS 



i' 
Qk 



■W R v-vec{AU s }}+0{A 2 } 
Im {t-^sls ' vec{A£/ s }} + O {A 2 } , where 

J 2 

n 



T 

Qk 



pI ■ (Ji ■ u s y 



T (Jl -U,' R1 

Pk 



Jl 

Fsls • F 



SLS 



(100) 



Xr 



a(p) 



(101) 



where a € 

rnNxl 



^Mxl 



is the array steering vector and s € 
contains the source symbols. Let Pj- — \\s\\? /N be 
the empirical source power. Furthermore, since we assume 
a ULA of isotropic elements, a(/i) is given by a(/j,) 



[l,e 



>e (Af-i)jM]. Note that || o-(m) || 2 = For 



notational convenience, we drop the explicit dependence of 
aon/i and write just a(/i) = a in the sequel. The selection 
matrices Ji and J 2 are then chosen as 



Ji = [J M -i M -ixi] J 2 = [0 M -ixi Im-i] (102) 

for maximum overlap, i.e., = M - 1. Since fllOll i is 

a rank-one matrix, we can directly relate the subspaces to the 
array steering vector and the source symbol matrix, namely 

TT a 1 

U s = u R = — — 



V 



M 
1 



(103) 
(104) 
(105) 

For the MSE expression from Theorem [3] we also require the 
quantity U n ■ J7„ , which resembles a projection matrix on th 
noise subspace. However, since the signal subspace is spanned 



NI2 ^P T -N 
M-N-Pt. 



\\a 



by a we can write U n ■ U n = Pr a = Im 

jj ■ a ■ a H . The MSE expression for 1-D Standard ESPRIT 
also include the eigenvectors of * which for the special case 
discussed here is scalar and given by \& = c m . Consequently, 
we have p k = q k = 1 for the eigenvectors. 

Combining these expressions and inserting into the special 
case of d33l for white noise, which is shown in equation 

d74l ), we have E 
o*/2-\\r T -W m 

(r) 

r k r = 

w mat 



{(A^) 2 } = a 2 /2 



W 



with 



Ji 



(J 2 /e^-J 1 ) 



Pri 



yjM ■ N ■ Pj ^/Pt ■ N. 



(106) 



(107) 



Note that W ma t is the Kronecker product of a 1 x N vector 
and an M x M matrix. Hence, r T • W ma t can be written as 



W, 



S < 



a , 



where 



which is the desired result d42l . Equation d43l follows 
from (f42T) by inserting the first order expansion for vcc { A?7 S } 
in terms of vec {TV} as shown in Appendix Icl □ 

Appendix F 
Proof of Theorem[6] 

This theorem consists of several parts which we address in 
separate subsections. 



- T 

a 



\Ji\l ■ N 
Ji 



Pl 



\/P T -N 



Ji ■ R- 



(108) 



(109) 



Therefore, the MSE can be expressed as E |(A/^[. r ' ) ) 2 | = 

cr 2 /2- ||s T <g> a T \\ 2 r which is equal to E |(A^ r) ) 2 | = er 2 /2- 
1 , 2 



|5 T I 



I ~T 

a 
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Since s T is a scaled version of s H and ||s H ||2 = N ■ we 
find that the first term in the MSE expression can conveniently 
be expressed as 



|S T I 



1 



P T -N 



1 



M ■ N ■ P T P T -N M ■ N -P T 



(HO) 



T 

Next, we proceed to simplify a further. To this end, we 
expand the pseudo-inverse of J x a using the rule x + = 
cc H / \ \x\\1 and multiply the brackets out. After straightforward 
algebraic manipulations we obtain 



-T 

a = 



<M 



-T 



-T 



where 

aj = a H • ■ J 2 /e J " and a J 

Since we have a H = [1, e~ m , e~ 2m , . . 
to see that 



= a H -jf -Ji. (Ill) 



-T 

a, 



-1 

a., 



[0,e- J ",e- 2j ", 



e -(M-2) m e -(M-l)j/il 



[1, e"^, e" 2 ^, . . . , e -( M - 2 )^, 0] and hence 



-1,0, 



O.e-^- 1 ^]. 

2 



a 



Consequently, we find 
result with (11 101 ) we finally 

e{(a m M) 2 } = | 



M 

L - (M-l)* ' 

have 



(112) 

2. Combining this 



l-TII 2 
1° h 



2 M-N-P T 
ol 1 



• 2 



M 



(M -1)2 



iV.p T (A/-1) 2 
which is the desired result. 



(113) 

(114) 
□ 



B. MSE for Unitary ESPRIT 

The second part of the theorem is to show that for a single 
source, the MSE for Unitary ESPRIT is the same as the MSE 
for Standard ESPRIT. Firstly, we expand Xq and find 



X 



(fba) 



[a ■ s T , Um ■ a* ■ s H ■ II jv] 



= a 



-,H 



n 



— T 

as 



(115) 



where we have used the fact that for our ULA we have II m 
a* = a - q-J^' 1 - 1 ) and we have defined s to be 



since Forward-Backward Averaging destroys the circular sym- 
metry of the noise, which leads to an additional term in the 
MSE expressions. Following the lines of the derivation for 1-D 
Standard ESPRIT we can show that 



r (fbar. w (fba) =i 



-T 



-T 

a , 



where s is given by 



-N -P T \/2 ■ N ■ P T 



(118) 



(119) 



and d is the same as in the derivation for 1-D Standard 
ESPRIT (cf. equation <fl09l ), 

According to Theorem |4] the MSE for Unitary ESPRIT can 
be computed as 



-^■(Z T -Z* -R C {Z T -U2MN-Z}), 



(120) 



(fba 



^mat > where we have inserted R^ a 



CT n ' J-2MN and = a 2 ■ TL2MN since we are consid- 

ering the special case of circularly symmetric white noise. 
Using (1118b and the fact that U 2 mn = n 2 jv ® Hm> this 
expression can be written into 



-II 2 11 — 11 2 
2 v M^lla H«-ll- 



S ■ U-2N ■ 8 -Ob ■ IljVf • d ) . (121) 



Since d is the same as in ( II 1 II ) we know that 1 1 ce. 1 1 ^ = pjzjp-- 



11— n 2 
Moreover, l|s|| 



— - follows directly from ( 11 19b . For 
the second term in ( 1121b we have 



s ■ Y1 2 n ■ s 
1 



1 



2MNP T 2NP T 
1 



2MNP T 2NP T 

1 . e w(Af-i) 
2MNP T 



oJA i(A/-l) 



(122) 



Similarly we can simplify d T -Hm-o, by using (II 1 lb and ( II 12b . 
We obtain 



d ■ Um ■ d 



2M 



(M-iy 



-(M-X)j/x 



(123) 



e - m {M-l) . jj n . s * 



(116) 



Note that s H s = s H 



s + s L ■ n N ■ n N ■ s* = 2 ■ s n ■ s. As 
for Standard ESPRIT, we relate <fTT3b to its SVD and obtain 



,(fba) 



.(fba) 



a 



.(fba) 



s 



>M 



\/2 ■ N ■ P T ' 



2 ■ M ■ N ■ P T . 



(117) 



An important consequence we can draw from (II 17b is that 
the column space u s remains unaffected from the forward- 
backward-averaging. Therefore we also have [/f 5a) = U n and 
hence U { ^ ] U^ ha)H = I M - ^ft- However, the equivalence 
of Standard ESPRIT and Unitary ESPRIT is still not obvious 



Combining the results from (1122b and (1123b into (1121b we 
finally obtain for the MSE 



1 



2M 



fi (— 

2 ' \2MNP T ' (M- l) 2 

1 ,,,C)i/r_Ti 2M 



2MNP T 
' 1 



1 



(M - l) 2 
1 



-(M-1) JM 



1 



2 \NP T (M-l) 2 N P T (M-l) 2 
^ 1 



N .p T (M-l) 2 ' 



(124) 
(125) 



which is equal to the result for 1-D Standard ESPRIT 
from ( II 14b and hence proves this part of the theorem. □ 
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C. Cramer-Rao Bound 

The third part of the theorem is to simplify the deterministic 
Cramer-Rao Bound (CRB) for the special case of a single 
source. To this end, a closed-form expression for the deter- 
ministic CRB for this setting is given by ll34l 

.2 , r n .rp.-l 



c = 



2-N 



■ Re- 



D Pr^ D 



(126) 



where i?s = jf ■ S • S H is the sample covariance matrix of the 
source symbols, D € C Mxd is the matrix of partial derivatives 
of the array steering vectors with respect to the parameters of 

interest, and Pr^ = I M - A ■ • Aj ■ A H . In the case 

d — 1, we have i?s = ||s|| 2 /N = P\ and the CRB expression 
simplifies into 

-l 



C 



2-N -Br 



Rc i d 



■M 



M 



d 



1 

2~p 
1 

2~p 



Rc <M • d 



M 



■dT-a-a H -d 



d n d 



M 



\d H -a\ 2 



(127) 
(128) 
(129) 



Since for a ULA the array steering vector can be expressed as 
a = [l e"* e 2 ^ . . . e^- 1 ^"] we have 



da r „ 
d = — = j ■ 2 • c 2 ^ 

dp L 



(M — 1) ■ e ( M " 1 )»] . 

(130) 



Consequently the terms d H • d and d H • a become 

M-l 

d H • d = ^2 m 2 = - ■ (M — 1) ■ M ■ (2M - 1) (131) 



m=0 



G 

M-l 
m=0 



-j.-.(M-l)-M (132) 



Using these expressions in ( 11291 ), we obtain 



C 



2-p 
1 

~2~p 
1 



- • (M — !)■ M ■ (2M - 1) 



1 

M 



12 



3- I- (M-l) -M 2 



(M - 1) ■ M ■ (M + 1) 
6 



p (M-l)-M-{M + l)' 
which is the desired result. 



(133) 

(134) 
□ 



Appendix G 
Proof of Theorem[7] 

As shown in Theorem [5] the MSE for SLS in the special 
case of circularly symmetric white noise can be expressed as 

2 

e{(A^ sls ) 2 } = ^ • ||W mat • r^sLsH 2 , , (135) 
where r k . SLS = r k r S - Ar k . SLS , 



Ar 



fe,SLS 



T {Jl-U S ) K 



SLS 



and Wr.u as well as -Fsls are defined in Theorem [5] For 
a single source, we have p k = q k = 1, \l> = ^ = e"*, and 
C/ s = aj\[M, and therefore Ar^ gLg = Ar gLg simplifies to 

(Ji • a) H 



Ar T 
^ r SLS 



' SLS ' -FsLS 



W 



R,U 



(136) 



Wr,u = (e^ • Ji) + [Ji ■ a ( Ji -ay -J 



SLS 



Jl 



Ji • a ( Ji -a) ■ Ji I — J 



, e 3 " Ji J 2 



M 



We can write IVr u as 

W R>V = ( J i ■ a ( J i ■ a) + - I M _i 



(Ji 



3 " ■ Ji) 
(137) 

Moreover, we need to simplify the term ^Fsls ■ F gLg ^ ■ ^ 
is easily verified that .Fsls ■ Fg Lg can l> e written as 
^sls ■ Fsls = diag ( J x • a) • G ■ diag ( J x ■ af (138) 

G = — • l(M-l)x(M-l) + 2 • Jiif-i - Ji ■ J% - J2 ■ J\ 

(139) 



Equation (1138b shows that the inverse of .Fsls 1 F SLS can be 
expressed as diag ( J\ ■ a) ■ G" 1 ■ diag (J\ ■ a) H . To proceed 
further, we require the following Lemma: 

Lemma 3. 77ze inverse of the matrix G defined in J 1 39b is 
given by the following expression 

TT ■ [{M -mi) -m 2 



i G ~ 1 ] i 

I J (mi 



-3 ■ 



m-i • ( M — mi ) -7712 -(iVf — TU'z) 



(m 1 ,m 2 ) 



= < 



A/ 2 + ll 

h ■ ( m i • ( M ~ TO 2) 



-3- 



mi ■ ( M— mi ) -m2 -(M— 1TI2) 
M' 2 + ll 



mi > 1712 

mi < mi 
(140) 



mi, ?ri2 = 1, 2, . . . , M — 1. 



Proof: To prove this Lemma it is sufficient to multiply 
G 1 in ( 1140b with G defined in ( 11391 ) and show that the result 
is an identity matrix. ■ 
Collecting our intermediate results from ( 1136b , ( 11381 ). and 
(1137b we have for Arg Lg 

(./ia) H h 
^ r SLS = — ■ diag (J i ■ a) ■ G 1 ■ diag ( Ji ■ a) 



'M ■ e-" 1 
_ 7 (M) 



M 



(Ji ■ a) + ■ (J 2 /e"« - Ji) 



T 

diag (a) H 



(141) 



where the scalar 7(M) and the row-vector are defined as 



7 (M) 



L lx(M-l) 



G 1 • 1(M-1 



)xl 



3d = lix(M-i) • G 1 -(J 2 -Ji)e 



nlxM 



(142) 
(143) 
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For 7(M) we can show via Lemma [3] 



(mi,m2 



(M- 1)M(M + 1) 
M 2 + 11 ' 



M-l M-l 

7(M) = E E 

mi — 1 m2-l 

Moreover, for the m-th element of the vector g^, which we 
denote as <?D,m we can show 

9B, m = A/2 6 +n • (2m - M - 1), m = 1,2, . . . , M 

(144) 

Collecting our intermediate results, we have shown that rg LS 
can be written as 



SLS 



LS 



1 T j- / \H 

— • g v ■ diag (a) 



(145) 



where t"ls has been taken from ( 1106b . The next step to 
computing the mean square error is to calculate the squared 
norm of the vector W^ at • tsls- The first few steps in 
computing this product are very similar to the LS case. 
Following ( 1109b we find that in the SLS case, the result is 
again equal to the product of the squared norm the same 



i.e. 



and 



SLS 



SLS 



■W r 





1 








-T 






"SLS 




2 





where 



l SLS 



= r 



SLS 



Pr 



(146) 



has been computed in ( 1145b . Applying similar 

T 



arguments as in ( 111 11 1. a SLS can be further simplified into 



~ T 



a 



SLS 



12 ' (M - 1)(M 2 
M 6 



11) 



-(M-l)jm 



M 
-MH 



M 2 + 11 

l),(-M + 3)e-*\ 



[-l,0,...,0,e 



,{M- 1) . e ~ {M - lhfi ] 



We conclude that the two vectors a SLS consists of have the 
same phase in each element and can hence be conveniently 
combined. When computing the squared norm of a SLS by 
summing the squared magnitude of all elements the phase 
terms cancel which also confirms the intuition the the result 
should be independent of the particular position [i. We obtain 



l SLS 



'M 

M-l 

E 

m=2 



-12 



(M- 
3G 



M(M 2 



1)(M 2 

(-M 



11 



=12- 



M 
M 4 



ll) 2 

12 



6- (-M + 1) 
+ M(M 2 + 11) 

2m- l) 2 

6 • (M - 1) 



(M 
2M 3 



1)(M 2 
24M 2 - 



-11) 
22M 



M(M 2 
23 



11), 
(147) 



(M 2 + 11) 2 (M- l) 2 
The mean square error is given by (cf. equation ( II 13l >) 
E{(A / i SLS ) 2 } = 4 
and ( 1147b we have 

E{(A^ SLS ) 2 } 



-T 
"SLS 



12 M 4 - 2M 3 + 24A'/ 2 - 22M + 23 



2 MNP T 



(M 2 + 11) 2 (M -If 



• 6 



M 4 - 2M 3 + 24A/ 2 - 22M + 23 



(148) 
□ 



N-P T M(M 2 + 11) 2 (M — l) 2 

which is the desired result. 

Appendix H 
Proof of Theorem[8] 

A. R-D Standard ESPRIT 

The proof for the R-D extension is in fact quite similar to 
the proof for the 1-D case provided in Section[F] In fact, ( 1 1 lb 
is still valid, the only difference being that a(/j) becomes 
a(fj,^') ® . . . ® a(iu,^) = a(fi). Therefore, the first steps 
of the derivation can still be performed in the very same way. 
We obtain the MSE for R-D Standard ESPRIT as 



(r) 



2 



2 



Silo 



H f (r)H 

a n J 1 

fir) 

J 1 a 



(r) 



Ji 



(r) 



Pr 



-(r) 

Since J 1 selects the M r 



1 out of AL elements in the r- 



th mode, we have 



a 



M ■ M r 



M ■ (M r - 1) 

- a H.jT )H 



aW®...a( fl ) and ft 



for £=1,2 and r = 1,2,...,/?, all "unaffected" modes can 
be factored out of ( 1151b and we have 



M- M r 



M ■ (M r - 1) 



a 



(i) 



where and a 2 are given by 



~(r) J 

a 2 



17 1 

(r) H t(0 H 
' ** 1 



(r) 



,.(«) 



and 



J 1 



Inserting (QJ0) Following the same rea soning as for ( fTT2b we find 



-4 V =[-l,0,...,0,e-^- 1 ^] 



(149) 



where s is the same as in the 1-D case (cf. equation ( 1109b ) 
and a} r > is given by 



(150) 



= jj- ■ (M r - 1). Moreover, 
multiplying ( 1150b out and using the fact that a satisfies the 
shift invariance equation in the r-th mode, we obtain 



(151) 



Since the the array steering vector a and the selection matrices 

— (r) 

J \ can be factored into Kronecker products according to a = 



(152) 



(153) 



(154) 
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Consequently, the desired norm 
be 

r-1 

^'(ni» 



is directly found to 



a< r > 



M ■ Ml I W II (n) 5 



2 A/ 2 • (M r -iy 



=2 • 



2- J] ||a<"> 

n— r+1 

(155) 



With the help of ( 11581 ) we find for the diagonal elements (r\ 

r 2 = r) 

2 ^>\ d *> = «.fi>A 



M r 



m=0 



Therefore, the MSE expression for R-D Standard ESPRIT is 
given by 

- 2 1 „ Mr 



M ■ (M r - 1) • (2M r - 1) (160) 



and similarly 



• 2 • 



2 M-N-Pt (M r - 1) 



f/ ■ a = xr ■ -jL' 



Ar.p T M-(Af r -l) 2: 
which proofs the first part of the theorem. 



(156) 
□ 



1 



B. R-D Unitary ESPRIT 

The second part of the theorem is to prove that in the R- 
D case the performance of R-D Unitary ESPRIT and R-D 
Standard ESPRIT are the same as long as a single source is 
present. However, for this part, no changes have to be made 
compared to Appendix IF-B I As it was shown there, Forward- 
Backward-Averaging only affects v s and has no effect on u s 
or U n . Applying the same steps here immediately proves this 
part of the theorem. 

C. Cramer-Rao Bound 

The third part is the simplification of the Cramer-Rao 
Bound. In the R-D case, the CRB is given by 

^2 



= -J ■ M ■ - ■ (M r - 1). 

Combining these two results we have for [J] ^ r r ^ 

[J] {r>r) = ~M-(M r -l)(M r + l). 



(161) 



(162) 



On the other hand, for the off-diagonal elements we obtain 

df ri) ■ d [r2) = -M ■ (M ri - 1) • M ri ■ (M r2 - 1) • Mr 2 
and therefore for [J]( ri ,r 2 )> r i ^ r 2 



[J] 



■Jin) 11 r(r 2 ) -(n) H h T(' r 2) H _ 

-a -a —a a - a a = 0. 



This shows that J is diagonal and real-valued. Consequently, 
the CRB becomes 



C = 



2 • N 
where Pr] 



•Re 



'■M 



A* -A 



A H and DW G 



C (r) = 



2 ■ N ■ P T 



Re{J} _1 =diag( 



C (R) 



12 



(r^Mx(d-R) con t ams th e partial derivatives of the array steering 
vectors a n with respect to [±r^ for n — 1,2, ...,d and 
r = 1, 2, . . . , R. For the special case of a single source, the 
CRB simplifies into (cf. Appendix IF-Cb 

C= -Re{jr\ 



2-N-P T M- (M r - 1) • (M r + 1) 
1 6 



(163) 
□ 



2 ■ N ■ P T 



J = ■ ( I M - - ■ a ■ a H ) • 



M 



(157) 



The columns of D {R) eC MxR are given by dt r) = e 



iMxl 



Using the fact that a = a 



(i) 



a {R) we obtain 



,(*) 
(158) 



where d 



gqW 
9^ 



iM r Xl 



J • [0,e^ M ,2 



e 2 ^ , . . . , (M — l)e( Af - 1 )^ ]. Therefore, the elements of 
the matrix J are given by 

[ J ](n,r 2 ) = d d ~M aa d 

(159) 



p M-(M r -l)-(M r + l)' 
which is the desired result. 

D. R-D Standard Tensor-ESPRIT 

The fourth part of the theorem is to show that the MSE of 
R-D Standard Tensor-ESPRIT is the same as the MSE for R- 
D Standard ESPRIT for d = 1. Since we have only shown the 
expressions for R-D Standard Tensor-ESPRIT in the special 
case R = 2, we will also assume this case here. 

Note that the MSE expression for Tensor-ESPRIT is in fact 
quite similar to the one for matrix-based ESPRIT with the 
only difference being that the matrix W ma t is replaced by 
the matrix Wtem cf- (f35T > and d36*l l. respectively. 

To simplify this expression for the special case d — 1, we 
express the unfoldings of Xq as 



X Q ] {1) = a w • [a (2) ® sj , [X a ] {2) = a (2) • 

_„ f-.(l)^„(2A T 

(3) 



(1) 



[Xo\,^ = s ■ | a v "' C$ a 
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Consequently, we can relate the necessary subspaces of the 
unfoldings of Xq to s and via 



a 



(i) 



,<2) 



U 



jM _ p,.-L r r[n] _ „ J_ 



E^ 1 = E 2 s] = Eg' = \IM ■ N ■ P T 



v\ = 



(a 



(2) 



s <g> a 



(164) 



a/m 2 - n- P t ' 

[b] a 



tfi* • vr r = u a ■ u» = 

Moreover, we have for T r 

T r = u [ * ] ■ u [ * ]H = Pr a(r) for r = 1, 2 and thus 
T l( g>T 2 =PP a( D ®Pr a( 2) = Pr a . (165) 

From (1164b and dl65l l it immediately follows that the first term 
in Wton cancels as it contains [Ti ® T 2 ] • V [ £ ] * ■ V [ £ ] . We 
also find *„„ = ■ e"^ M /M r for r = 1,2 and 
m = 1,2, ...,M r . To simplify the remaining two terms in 
Wton we first simplify some of their components. Using the 
identity u s = a/\M and the explicit expression for t r<m it 
is easy to show that 



(uj ® I M ) • Ti = — L fa (2)T ® a (1) ® J M „ 
(uj (g) J M ) • T 2 = — L fa (1)T <g) I Ml ® a (2) 



(166) 
(167) 



Moreover, using the relations from ( 1164b we can rewrite 



as 



i 



MNP T 



,( 2 ) 



2 ^2 

=^^- • f a< 2 )* • (« g) a«) H ) ® P^ (2) . 
MNP T V v ; J 

Combining these intermediate result, the third term in Wt cn 
can be expressed as 

K T ®i M ).iv((i/r srvi 



Af 2 



s H «ia (1)H 



,.(!) 



3 {u [ ?u™ 

(168) 



MTVPtVM 

With similar arguments, the second term in Wt cn can be 
simplified into 



Wt 



Using (11681 1 and (1169l l in (f36]l, we obtain 
1 



NP T VM 



(170) 



Comparing (1170b matrix-based counterpart in ( 1150b we find 
that for a single source, W mat and Wtcn are in fact quite 
similar, the only difference being that Pr^ is replaced by 
Pr„ ( i) <2> B" a C2) + R" a (i) <£> R"^ ( 2) ■ Therefore, to show the R- 
D Standard ESPRIT and R-D Standard Tensor-ESPRIT have 
the same MSE for d = 1, it is sufficient to show that the 
corresponding terms ar- r ' are the same, i.e., that 

a« • Pr^ = a« • (Pr^ (1) ® Pr a(2) + Pr a(1) ® Pr^ (2 



(171) 

j[ r) ) for r = 1,2. 
Note that • Pr^ was shown to be equal to (cf. equa- 

T 



where a« = a^j^ (j^ / e ^ ( " 



tion (1152b ) lai 



=(i) 



for 



1 and a' 1 ' 



(r) H 



(a< 2) - af) for r = 2, where a< r)T = aM H • jj 

J 2 r) /c J7t<r> and a 2 r)T = a< r ) H • jf' )H • jj r) . Expanding the 
corresponding right-hand side of ( 1171b we have for r = 1 



" (r) ' few ® + ® ^at 2 ) 



= (a( 1 ) T -a< 2 ) T )®a( 2 ) H , 

where we have used the fact that a = a^ 1 ) <g>a,( 2 ) from its defi- 
nition, a ( 2 ) -Pr a (2) = a^ 2 ^ and ■Pr a ( 2 ) = 0i x m 2 since 
Pr a( 2) and Pr (2 ) are projectors onto a' 2 ) and its orthogonal 
complement, and the identity J^'/e 1 '^ — J^a^ = 

0(M 2 -i)xi since a' 1 ) satisfies the shift invariance equation for 
r = 1. This shows that the left-hand side and the right-hand 
side of (1171b are equal for r = 1. The proof for r = 2 proceeds 
in an analogous fashion. Consequently, we have shown that for 
d=l 



r {r)T ■ W mat = r« T • W tcn , for r = 1, 2 



(172) 



Af 2 x(MfJV) 



and hence the MSE for 2-D Standard ESPRIT and 2-D 
Standard Tensor-ESPRIT are in fact equal. □ 

E. R-D Unitary Tensor-ESPRIT 

The fifth and final part of the theorem is to show that the 
MSE for R-D Unitary ESPRIT is again equal to the MSE 
for R-D Standard ESPRIT in case of a single source. Again, 
there is no need to derive this in full detail. As it was shown in 



MNPtVM 
Mi 



,(2) H 



RaU> 



,12) 



,(2: 



,(2) H 



K 



J\/ 2 x(Mi-JV) 

(169) 



MNP T VM 

where the last step is a special case of Property ([H} for 
commutation matrices. 



Appendix IF-B I Forward-Backward- Averaging has no effect on 
tt s or U n but only affects v s and V n . This carries over to the 
tensor case where only the quantities involving the symbols are 
affected. However, since the "symbol part" and the "array part" 
can always be factorized (cf. equation ( 11701 )). the arguments 
from Appendix IF-BI can still be applied to prove this part of 
the theorem. □ 



24 



References 

[1] T. W. Anderson, "Asymptotic theory for principal component analysis", 

Ann. Math. Statist., vol. 34, pp. 122-148, 1963. 
[2] D. R. Brillinger, Time Series: Data Analysis and Theory, Holt, Rhinehart 

and Winston, New York, 1975. 
[3] L. De Lathauwer, "First-order perturbation analysis of the best rank- 

(r i , T2 , T3 ) approximation in multilinear algebra", Journal of Chemo- 

metrics, vol. 18, no. 1, pp. 2—11, 2004. 
[4] L. De Lathauwer, B. De Moor, and J. Vandewalle, "A multilinear 

singular value decomposition", SIAM J. Matrix Anal. AppL, vol. 21, 

no. 4, pp. 1253-1278, 2000. 
[5] J. R Delmas and Y. Meurisse, "On the second-order statistics of the 

EVD of sample covariance matrices - application to the detection of 

noncircular or/and NonGaussian components", IEEE Transactions on 

Signal Processing, vol. 59, no. 8, pp. 4017-4023, Aug. 2011. 
[6] A. Eriksson, P, Stoica, and T. Soderstrom, "Second-order properties of 

MUSIC and ESPRIT estimates of sinusoidal frequencies in high SNR 

scenarios", IEE Proceedings-F , vol. 140, no. 4, Aug. 1993. 
[7] B. Friedlander, "A sensitivity analysis of MUSIC algorithm", IEEE 

Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP- 

38, pp. 1740-1751, 1990. 
[8] M. Haardt, "Structured least squares to improve the performance of 

ESPRIT-type algorithms", IEEE Transactions on Signal Processing, 

vol. 45, pp. 792-799, Mar. 1997. 
[9] M. Haardt and J. A. Nossek, "Unitary ESPRIT: How to obtain increased 

estimation accuracy with a reduced computational burden", IEEE 

Transactions on Signal Processing, vol. 43, no. 5, pp. 1232-1242, May 

1995. 

[10] M. Haardt and J. A. Nossek, "Simultaneous Schur decomposition 
of several non-symmetric matrices to achieve automatic pairing in 
multidimensional harmonic retrieval problems", IEEE Transactions on 
Signal Processing, vol. 46, no. 1, pp. 161-169, Jan. 1998. 

[11] M. Haardt, F. Roemer, and G. Del Galdo, "Higher-order SVD based 
subspace estimation to improve the parameter estimation accuracy in 
multi-dimensional harmonic retrieval problems", IEEE Trans. Sig. Proc. , 
vol. 56, pp. 3198 - 3213, July 2008. 

[12] M. Haardt, R. S. Thoma, and A. Richter, "Multidimensional high- 
resolution parameter estimation with applications to channel sounding", 
in High-Resolution and Robust Signal Processing, Y. Hua, A. Gershman, 
and Q. Chen, Eds., pp. 255-338. Marcel Dekker, New York, NY, 2004, 
Chapter 5. 

[13] Y. Hua and T. K. Sarkar, "On SVD for estimating generalized 
eigenvalues of singular matrix pencil in noise", IEEE Transactions on 
Signal Processing, vol. 39, pp. 892-900, Apr. 1991. 

[14] M. Kaveh and A. J. Barabell, "The statistical performance of the MUSIC 
and the minimum-norm algorithms in resolving plane waves in noise", 
IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 
ASSP-34, pp. 331-341, Apr. 1986. 

[15] H. Krim, P. Forster, and J. G. Proakis, "Operator approach to perfor- 
mance analysis of root-MUSIC and root-min-norm", IEEE Transactions 
on Signal Processing, vol. 40, no. 7, pp. 1687-1696, July 1992. 

[16] P. Lancaster and M. Tismentsky, The Theory of Matrices, New York 
Acadmic Press, 2nd edition, 1978. 

[17] F. Li, H. Liu, and R. J. Vaccaro, "Performance analysis for DOA 
estimation algorithms: Unification, simplifications, and observations", 
IEEE Transactions on Aerospace and Electronic Systems, vol. 29, no. 
4, pp. 1170-1184, Oct. 1993. 

[18] F. Li and R. J. Vaccaro, "Performance degradation of DOA estimators 
due to unknown noise fields", IEEE Transactions on Signal Processing, 
vol. 40, no. 3, pp. 686-690, Mar. 1992. 

[19] J. Liu, X. Liu, and X. Ma, "First-order perturbation analysis of singular 
vectors in singular value decomposition", IEEE Transactions on Signal 
Processing, vol. 56, no. 7, pp. 3044-3049, July 2008. 

[20] X. Liu, N. D. Sidiropoulos, and A. Swami, "Blind high-resolution 
localization and tracking of multiple frequency hopped signals", IEEE 
Transactions on Signal Processing, vol. 50, no. 4, pp. 889-901, Apr. 
2002. 

[21] J. R. Magnus and H. Neudecker, Matrix differential calculus with 
applications in statistics and econometrics, John Wiley and Sons, 1995. 



[22] C. P. Mathews, M. Haardt, and M. D. Zoltowski, "Performance analysis 
of closed-form, ESPRIT based 2-D angle estimator for rectangular 
arrays", IEEE Signal Processing Letters, vol. 3, pp. 124-126, Apr. 
1996. 

[23] C. P. Mathews and M. D. Zoltowski, "Performance analysis of the 
UCA-ESPRIT algorithm for circular ring arrays", IEEE Transactions 
on Signal Processing, vol. 42, no. 9, Sept. 1994. 

[24] D. Nion and N. D. Sidiropoulos, "Tensor algebra and multidimensional 
harmonic retrieval in signal processing for MIMO radar", IEEE 
Transactions on Signal Processing, vol. 58, no. 11, pp. 5693-5705, Nov. 
2010. 

[25] B. Ottersten, M. Viberg, and T. Kailath, "Performance analysis of the 
Total Least Squares ESPRIT algorithm", IEEE Transactions on Signal 
Processing, vol. 39, no. 5, pp. 1122-1135, May 1991. 

[26] S. U. Pillai and B. H. Kwon, "Performance analysis of MUSIC- 
Type high resolution estimators for direction finding in correlated and 
coherent scenes", IEEE Transactions on Acoustics, Speech, and Signal 
Processing, vol. AASP-37, no. 8, Aug. 1989. 

[27] B. Porat and B. Friedlander, "Analysis of the asymptotic relative 
efficiency of the MUSIC algorithm", IEEE Transactions on Acoustics, 
Speech, and Signal Processing, vol. 36, no. 4, pp. 532-544, Apr. 1988. 

[28] B. D. Rao and K. V. S. Hari, "Performance analysis of ESPRIT and 
TAM in determining the direction of arrival of plane waves in noise", 
IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 
AASP-37, no. 12, Dec. 1989. 

[29] B. D. Rao and K. V. S. Hari, "Performance analysis of Root-MUSIC", 
IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 
37, no. 12, pp. 1939-1948, Dec. 1989. 

[30] F. Roemer, H. Becker, M. Haardt, and M. Weis, "Analytical performance 
evaluation for HOSVD-based parameter estimation schemes", in Proc. 
of the IEEE Int. Workshop on Comp. Adv. in Multi-Sensor Adaptive 
Proc. (CAMSAP 2009), Aruba, Dutch Antilles, Dec. 2009. 

[31] R. Roy, A. Paulraj, and T. Kailath, "ESPRIT - a subspace rotation 
approach to estimation of parameters of cisoids in noise", IEEE 
Transactions on Acoustics, Speech, and Signal Processing, vol. AASP- 
34, pp. 1340-1342, Oct. 1986. 

[32] R. O. Schmidt, "Multiple emitter location and signal parameter estima- 
tion", IEEE Transactions on Antennas and Propagation, vol. AP-34, 
no. 3, pp. 243-258, Mar. 1986. 

[33] S. T. Smith, "Statistical resolution limits and the complexified cramer- 
rao bound", IEEE Transactions on Signal Processing, vol. 53, no. 5, 
pp. 1597-1609, May 2005. 

[34] P. Stoica and A. Nehorai, "MUSIC, maximum likelihood, and Cramer- 
Rao bound", IEEE Transactions on Acoustics, Speech, and Signal 
Processing, vol. 37, no. 5, pp. 720-741, May 1989. 

[35] P. Stoica and T. Soderstrom, "Statistical analysis of MUSIC and sub- 
space rotation estimates of sinusoidal frequencies", IEEE Transactions 
on Signal Processing, vol. 39, no. 8, Aug. 1991. 

[36] A. L. Swindlehurst, B. Ottersten, R. Roy, and T. Kailath, "Multiple 
invariance ESPRIT", IEEE Transactions on Signal Processing, vol. 40, 
no. 4, pp. 867-881, Apr. 1992. 

[37] A. Thakre, M. Haardt, and K. Giridhar, "Single snapshot r-d unitary 
esprit using an augmentation of the tensor order", in Prof. IEEE 
Int. Workshop on Computational Advances in Multi-Sensor Adaptive 
Processing (CAMSAP 2009), Dec. 2009. 

[38] A. Thakre, M. Haardt, and K. Giridhar, "Single snapshot spatial smooth- 
ing with improved effective array aperture", IEEE Signal Processing 
Letters, vol. 16, pp. 505-509, June 2009. 

[39] A. Thakre, M. Haardt, F. Roemer, and K. Giridhar, "Tensor-Based spatial 
smoothing (TB-SS) using multiple snapshots", IEEE Trans. Sig. Proc, 
Jan. 2010. 

[40] Z. Xu, "Perturbation analysis for subspace decomposition with appli- 
cations in subspace-based algorithms", IEEE Transactions on Signal 
Processing, vol. 50, no. 11, pp. 2820-2830, Nov. 2002. 

[41] M. D. Zoltowski, G. M. Kautz, and C. P. Mathews, "Performance 
analysis of eigenstructure based DOA estimators employing conjugate 
symmetric beamformers", in IEEE Sixth SP Workshop on Statistical 
Signal and Array Processing, Victoria, BC, Canada, Oct. 1992. 



